多層パーセプトロン (Multilayer perceptron, MLP)を scipy.optimize.curve_fit で実装する.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "scipy.optimize.curve_fit に関する記事を次々書いているところですが、\n",
    "\n",
    "* <a href=\"https://qiita.com/maskot1977/items/d16a385e23ed07e28037\">少ない観測値を補間してから、正規分布の線形和で近似する</a>\n",
    "* <a href=\"https://qiita.com/maskot1977/items/e4f5f71200180865986e\">カーブフィッティング手法 scipy.optimize.curve_fit の使い方を理解する</a>\n",
    "* <a href=\"https://qiita.com/maskot1977/items/0f8a1ac961845735ce06\">ロジスティック回帰を scipy.optimize.curve_fit で実装する</a> \n",
    "\n",
    "に続いて、今回は「多層パーセプトロン (Multilayer perceptron, MLP)を scipy.optimize.curve_fit で実装する」に挑戦します。\n",
    "\n",
    "# 多層パーセプトロン (Multilayer perceptron, MLP) とは\n",
    "\n",
    "過去に \n",
    "\n",
    "* <a href=\"https://qiita.com/maskot1977/items/5871e0b1b6a44f17c327\">多層パーセプトロン (Multilayer perceptron, MLP)をExcelで理解する</a>\n",
    "* <a href=\"https://qiita.com/maskot1977/items/d0253e1eab1ff1315dff\">多層パーセプトロン (Multilayer perceptron, MLP)をPythonで理解する</a>\n",
    "\n",
    "を書いたので、そちらをご覧ください。前者の方が内容は平易だと思います。\n",
    "\n",
    "# 「あやめのデータ」を読み込む\n",
    "\n",
    "「あやめのデータ」については <a href=\"https://qiita.com/maskot1977/items/453eddf5aa154c04668c\">実習用データ</a>　をご参照ください。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ウェブ上のリソースを指定する\n",
    "url = 'https://raw.githubusercontent.com/maskot1977/ipython_notebook/master/toydata/iris.txt'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# URL によるリソースへのアクセスを提供するライブラリをインポートする。\n",
    "import urllib.request\n",
    "urllib.request.urlretrieve(url, 'iris.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "df = pd.read_csv('iris.txt', delimiter=\"\\t\", index_col=0)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/kot/miniconda3/envs/py3new/lib/python3.6/site-packages/ipykernel_launcher.py:1: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    }
   ],
   "source": [
    "X = df.iloc[:, :4].T.as_matrix()\n",
    "Y = np.array(df.iloc[:, 4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 多層パーセプトロンの実装（出力層のノードが３つ）\n",
    "\n",
    "<a href=\"https://qiita.com/maskot1977/items/5871e0b1b6a44f17c327\">多層パーセプトロン (Multilayer perceptron, MLP)をExcelで理解する</a> と同じく、入力層Xに4つのノード、隠れ層Hには3つのノード、出力層Oに3つのノードを配置したMLPを実装しようと思います。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mlp(X, *params):\n",
    "    h01, h02, h03, h04, h0b = params[0], params[1], params[2], params[3], params[4]\n",
    "    h11, h12, h13, h14, h1b = params[5], params[6], params[7], params[8], params[9]\n",
    "    h21, h22, h23, h24, h2b = params[10], params[11], params[12], params[13], params[14]\n",
    "    \n",
    "    o01, o02, o03, o0b = params[15], params[16], params[17], params[18]\n",
    "    o11, o12, o13, o1b = params[19], params[20], params[21], params[22]\n",
    "    o21, o22, o23, o2b = params[23], params[24], params[25], params[26]\n",
    "    \n",
    "    h0 = 1. / (1. + np.exp(-(h01 * X[0] + h02 * X[1] + h03 * X[2] + h04 * X[3] + h0b)))\n",
    "    h1 = 1. / (1. + np.exp(-(h11 * X[0] + h12 * X[1] + h13 * X[2] + h14 * X[3] + h1b)))\n",
    "    h2 = 1. / (1. + np.exp(-(h21 * X[0] + h22 * X[1] + h23 * X[2] + h24 * X[3] + h2b)))\n",
    "    \n",
    "    o0 = 1. / (1. + np.exp(-(o01 * h0 + o02 * h1 + o03 * h2 + o0b)))\n",
    "    o1 = 1. / (1. + np.exp(-(o11 * h0 + o12 * h1 + o13 * h2 + o1b)))\n",
    "    o2 = 1. / (1. + np.exp(-(o21 * h0 + o22 * h1 + o23 * h2 + o2b)))\n",
    "    \n",
    "    return np.array([o0, o1, o2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... と、ここで手が止まった。「出力層のノードが３つ」だと...？\n",
    "\n",
    "今まで scipy.optimize.curve_fit で取り扱った関数のアウトプットは、スカラーだった。ベクトルをアウトプットする関数を scipy.optimize.curve_fit で取り扱えるのか？\n",
    "\n",
    "## 出力層に合わせて、Yも one-hot-vector にする\n",
    "\n",
    "そうやって scipy.optimize.curve_fit が取り扱ってくれればいいんだが..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/kot/miniconda3/envs/py3new/lib/python3.6/site-packages/ipykernel_launcher.py:1: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    }
   ],
   "source": [
    "Y = pd.get_dummies(pd.DataFrame([\"y=\" + str(y) for y in df.iloc[:, 4].as_matrix()])).as_matrix()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "Improper input: N=27 must not exceed M=3",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-85-da29f74a384b>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moptimize\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mcurve_fit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mpopt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpcov\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcurve_fit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmlp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mY\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp0\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# poptは最適推定値、pcovは共分散\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mpopt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/py3new/lib/python3.6/site-packages/scipy/optimize/minpack.py\u001b[0m in \u001b[0;36mcurve_fit\u001b[0;34m(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)\u001b[0m\n\u001b[1;32m    749\u001b[0m         \u001b[0;31m# Remove full_output from kwargs, otherwise we're passing it in twice.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    750\u001b[0m         \u001b[0mreturn_full\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'full_output'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 751\u001b[0;31m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mleastsq\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mDfun\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mjac\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    752\u001b[0m         \u001b[0mpopt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpcov\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minfodict\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0merrmsg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mier\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    753\u001b[0m         \u001b[0mcost\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minfodict\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'fvec'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m**\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/py3new/lib/python3.6/site-packages/scipy/optimize/minpack.py\u001b[0m in \u001b[0;36mleastsq\u001b[0;34m(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)\u001b[0m\n\u001b[1;32m    384\u001b[0m     \u001b[0mm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    385\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mn\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 386\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Improper input: N=%s must not exceed M=%s'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    387\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mepsfcn\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    388\u001b[0m         \u001b[0mepsfcn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfinfo\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meps\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: Improper input: N=27 must not exceed M=3"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit  \n",
    "popt, pcov = curve_fit(mlp,X,Y.T, p0=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]) # poptは最適推定値、pcovは共分散\n",
    "popt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ダメでした...(´・ω・`)\n",
    "\n",
    "# 多層パーセプトロンの実装（出力層のノードが１つ）\n",
    "\n",
    "出力層のノードが１つの場合は、今までと同じ。これまでの scipy.optimize.curve_fit に関する記事が理解できていれば、書けますよね！？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mlp(X, *params):\n",
    "    h01, h02, h03, h04, h0b = params[0], params[1], params[2], params[3], params[4]\n",
    "    h11, h12, h13, h14, h1b = params[5], params[6], params[7], params[8], params[9]\n",
    "    h21, h22, h23, h24, h2b = params[10], params[11], params[12], params[13], params[14]\n",
    "    \n",
    "    o01, o02, o03, o0b = params[15], params[16], params[17], params[18]\n",
    "    \n",
    "    h0 = 1. / (1. + np.exp(-(h01 * X[0] + h02 * X[1] + h03 * X[2] + h04 * X[3] + h0b)))\n",
    "    h1 = 1. / (1. + np.exp(-(h11 * X[0] + h12 * X[1] + h13 * X[2] + h14 * X[3] + h1b)))\n",
    "    h2 = 1. / (1. + np.exp(-(h21 * X[0] + h22 * X[1] + h23 * X[2] + h24 * X[3] + h2b)))\n",
    "    \n",
    "    o0 = 1. / (1. + np.exp(-(o01 * h0 + o02 * h1 + o03 * h2 + o0b)))\n",
    "    \n",
    "    return o0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "データも取り直し。<a href=\"https://qiita.com/maskot1977/items/0f8a1ac961845735ce06\">ロジスティック回帰を scipy.optimize.curve_fit で実装する</a> と同じものです。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "X1 = [4.7, 4.5, 4.9, 4.0, 4.6, 4.5, 4.7, 3.3, 4.6, 3.9, \n",
    "      3.5, 4.2, 4.0, 4.7, 3.6, 4.4, 4.5, 4.1, 4.5, 3.9, \n",
    "      4.8, 4.0, 4.9, 4.7, 4.3, 4.4, 4.8, 5.0, 4.5, 3.5, \n",
    "      3.8, 3.7, 3.9, 5.1, 4.5, 4.5, 4.7, 4.4, 4.1, 4.0, \n",
    "      4.4, 4.6, 4.0, 3.3, 4.2, 4.2, 4.2, 4.3, 3.0, 4.1, \n",
    "      6.0, 5.1, 5.9, 5.6, 5.8, 6.6, 4.5, 6.3, 5.8, 6.1, \n",
    "      5.1, 5.3, 5.5, 5.0, 5.1, 5.3, 5.5, 6.7, 6.9, 5.0, \n",
    "      5.7, 4.9, 6.7, 4.9, 5.7, 6.0, 4.8, 4.9, 5.6, 5.8, \n",
    "      6.1, 6.4, 5.6, 5.1, 5.6, 6.1, 5.6, 5.5, 4.8, 5.4, \n",
    "      5.6, 5.1, 5.1, 5.9, 5.7, 5.2, 5.0, 5.2, 5.4, 5.1]\n",
    "\n",
    "X2 = [1.4, 1.5, 1.5, 1.3, 1.5, 1.3, 1.6, 1.0, 1.3, 1.4, \n",
    "      1.0, 1.5, 1.0, 1.4, 1.3, 1.4, 1.5, 1.0, 1.5, 1.1, \n",
    "      1.8, 1.3, 1.5, 1.2, 1.3, 1.4, 1.4, 1.7, 1.5, 1.0, \n",
    "      1.1, 1.0, 1.2, 1.6, 1.5, 1.6, 1.5, 1.3, 1.3, 1.3, \n",
    "      1.2, 1.4, 1.2, 1.0, 1.3, 1.2, 1.3, 1.3, 1.1, 1.3, \n",
    "      2.5, 1.9, 2.1, 1.8, 2.2, 2.1, 1.7, 1.8, 1.8, 2.5, \n",
    "      2.0, 1.9, 2.1, 2.0, 2.4, 2.3, 1.8, 2.2, 2.3, 1.5, \n",
    "      2.3, 2.0, 2.0, 1.8, 2.1, 1.8, 1.8, 1.8, 2.1, 1.6, \n",
    "      1.9, 2.0, 2.2, 1.5, 1.4, 2.3, 2.4, 1.8, 1.8, 2.1, \n",
    "      2.4, 2.3, 1.9, 2.3, 2.5, 2.3, 1.9, 2.0, 2.3, 1.8]\n",
    "\n",
    "X3 = [7.0, 6.4, 6.9, 5.5, 6.5, 5.7, 6.3, 4.9, 6.6, 5.2, \n",
    "      5.0, 5.9, 6.0, 6.1, 5.6, 6.7, 5.6, 5.8, 6.2, 5.6, \n",
    "      5.9, 6.1, 6.3, 6.1, 6.4, 6.6, 6.8, 6.7, 6.0, 5.7, \n",
    "      5.5, 5.5, 5.8, 6.0, 5.4, 6.0, 6.7, 6.3, 5.6, 5.5, \n",
    "      5.5, 6.1, 5.8, 5.0, 5.6, 5.7, 5.7, 6.2, 5.1, 5.7, \n",
    "      6.3, 5.8, 7.1, 6.3, 6.5, 7.6, 4.9, 7.3, 6.7, 7.2, \n",
    "      6.5, 6.4, 6.8, 5.7, 5.8, 6.4, 6.5, 7.7, 7.7, 6.0, \n",
    "      6.9, 5.6, 7.7, 6.3, 6.7, 7.2, 6.2, 6.1, 6.4, 7.2, \n",
    "      7.4, 7.9, 6.4, 6.3, 6.1, 7.7, 6.3, 6.4, 6.0, 6.9, \n",
    "      6.7, 6.9, 5.8, 6.8, 6.7, 6.7, 6.3, 6.5, 6.2, 5.9]\n",
    "\n",
    "X4 = [3.2, 3.2, 3.1, 2.3, 2.8, 2.8, 3.3, 2.4, 2.9, 2.7, \n",
    "      2.0, 3.0, 2.2, 2.9, 2.9, 3.1, 3.0, 2.7, 2.2, 2.5, \n",
    "      3.2, 2.8, 2.5, 2.8, 2.9, 3.0, 2.8, 3.0, 2.9, 2.6, \n",
    "      2.4, 2.4, 2.7, 2.7, 3.0, 3.4, 3.1, 2.3, 3.0, 2.5, \n",
    "      2.6, 3.0, 2.6, 2.3, 2.7, 3.0, 2.9, 2.9, 2.5, 2.8, \n",
    "      3.3, 2.7, 3.0, 2.9, 3.0, 3.0, 2.5, 2.9, 2.5, 3.6, \n",
    "      3.2, 2.7, 3.0, 2.5, 2.8, 3.2, 3.0, 3.8, 2.6, 2.2, \n",
    "      3.2, 2.8, 2.8, 2.7, 3.3, 3.2, 2.8, 3.0, 2.8, 3.0, \n",
    "      2.8, 3.8, 2.8, 2.8, 2.6, 3.0, 3.4, 3.1, 3.0, 3.1, \n",
    "      3.1, 3.1, 2.7, 3.2, 3.3, 3.0, 2.5, 3.0, 3.4, 3.0]\n",
    "\n",
    "X = np.array([X1, X2, X3, X4])\n",
    "\n",
    "Y = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \n",
    "     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \n",
    "     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \n",
    "     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \n",
    "     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \n",
    "     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, \n",
    "     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, \n",
    "     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, \n",
    "     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, \n",
    "     1, 1, 1, 1, 1, 1, 1, 1, 1, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "結果表示の関数も <a href=\"https://qiita.com/maskot1977/items/0f8a1ac961845735ce06\">ロジスティック回帰を scipy.optimize.curve_fit で実装する</a> と同じ。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "def confusion_table(y_pred, y_true):\n",
    "    true_positives = [] # TP\n",
    "    false_positives = [] # FP\n",
    "    false_negatives = [] # FN\n",
    "    true_negatives = [] # TN\n",
    "    for y1, y2 in zip(y_pred, y_true):\n",
    "        if y1 >= 0.5:\n",
    "            if y2 >= 0.5:\n",
    "                true_positives.append(y1)\n",
    "            else:\n",
    "                false_positives.append(y1)\n",
    "        else:\n",
    "            if y2 >= 0.5:\n",
    "                false_negatives.append(y1)\n",
    "            else:\n",
    "                true_negatives.append(y1)\n",
    "    return (true_positives, false_positives, false_negatives, true_negatives)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "def show_result(TP, FP, FN, TN):\n",
    "    print(\"Accuracy: \", len(TP + TN) / len(TP + FP + FN + TN))\n",
    "    plt.figure(figsize=(12,4))\n",
    "    plt.hist([TP, FP, FN, TN], label=['TP', 'FP', 'FN', 'TN'], color=['blue', 'green', 'orange', 'red'])\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 計算実行"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.57166185e+01,  1.58250335e+01, -3.52251392e+00, -7.49692230e+00,\n",
       "       -5.89003641e+01,  3.80089003e-02,  2.37483939e-01,  2.87007844e-01,\n",
       "       -9.48869044e-01,  1.39257951e+00, -1.28133971e-01, -5.79757478e-01,\n",
       "        3.16485573e-01, -1.03528732e+00,  1.34439583e+00,  1.12009673e+02,\n",
       "       -5.14928444e+01, -5.51022639e+01, -3.22706767e+01])"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit  \n",
    "popt, pcov = curve_fit(mlp,X,Y, p0=[1,0,0,0,0, 0,0,0,0,0,\n",
    "                                     0,0,0,0,0, 0,0,0,0]) # poptは最適推定値、pcovは共分散\n",
    "popt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy:  0.99\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/kot/miniconda3/envs/py3new/lib/python3.6/site-packages/matplotlib/font_manager.py:1241: UserWarning: findfont: Font family ['IPAexGothic'] not found. Falling back to DejaVu Sans.\n",
      "  (prop.get_family(), self.defaultFamily[fontext]))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsMAAAD8CAYAAACSP6kTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAFWdJREFUeJzt3X+QXWd5H/Dvg+RYxpJlLFPhsWilGTsZfrheoR2GDGO6ayeUUGzhKeMxbYLpaKo/CgEau42p21q0wEBL44SZtI0bM6gdiEIZMBo3LvU4qEAnppGKcAANRaUiFTUxKMbxGpsg5+0fe+1KWPbe1d57V+v385nZ2XvOfe97Hs2jXX117nvPqdZaAACgR89b7gIAAGC5CMMAAHRLGAYAoFvCMAAA3RKGAQDoljAMAEC3hGEAALolDAMA0C1hGACAbq2e5MEuvPDCtnnz5rHN/+ijj+bcc88d2/ycWfS7P3reF/3ui373Zdz9PnDgwPdbay8cZuxEw/DmzZuzf//+sc2/b9++zMzMjG1+ziz63R8974t+90W/+zLuflfVt4cda5kEAADdEoYBAOiWMAwAQLcmumYYAIAzx49//OMcPXo0jz/++ESPu379+hw6dGjJ86xZsyabNm3KWWedddpzCMMAAJ06evRo1q1bl82bN6eqJnbcRx55JOvWrVvSHK21HDt2LEePHs2WLVtOe56hlklU1ZGq+qOqOlhV+wf7Lqiqe6rqm4PvLzjtKgAAmLjHH388GzZsmGgQHpWqyoYNG5Z8Vnsxa4ZnW2tTrbXpwfbNSe5trV2a5N7BNgAAK8hKDMJPGkXtS/kA3fYkuwePdyd545KrAQCACRp2zXBL8l+qqiX5rdba7Uk2ttYeGDz/3SQbx1EgAACTMeqTxK09+/PHjh3LVVddlST57ne/m1WrVuWFL5y/cdxXvvKVXH755Tl+/Hhe8pKXZPfu3Xn+858/2gKTVFuoyiRVdXFr7TtV9ZeS3JPkl5Psba2df8KYh1prT1s3XFU7k+xMko0bN27bs2fPyIr/SXNzc1m7du3Y5ufMot/90fO+6Hdf9PtkBw4s/jXbti3+NevXr88ll1zy1PZ55y3tQ20/6c/+7JFT7n/iiSeyatWqk/a9//3vz9q1a/OOd7wjSXLRRRflgQfmz7vu2LEjW7duzdvf/vanzXX48OE8/PDDJ+2bnZ09cMLS3mc11Jnh1tp3Bt8frKpPJ3llkj+pqotaaw9U1UVJHnyG196e5PYkmZ6ebuO89Z5bOfZFv/uj533R777o98lmZxf/miHObz7NoUOHlnxVh2fzTHOf6moSZ599ds4+++yT9j/5+Morr8z9999/yvnWrFmTrVu3nnaNC64Zrqpzq2rdk4+TvDbJV5PsTXLDYNgNST5z2lUAAMApHD9+PHfffXcuu+yyscw/zJnhjUk+Pfi03uokH2+t/eeq+sMkn6iqHUm+neS6sVQIAEB3HnvssUxNTSVJrrjiiuzYsWMsx1kwDLfWvpXk8lPsP5bkqnEU1Y3FrlI/nfc/AABWoHPOOScHDx4c+3GWcmk1AABY0dyOGQCAJH2+CS0MAwCw7Hbt2nXS9tzc3ESOa5kEAADdEoYBAOiWMAwAQLeEYQAAuiUMAwDQLWEYAIBuubQaAABJknrPIu+Ou4B268IXLl61alUuu+yyp7bvvPPOHDlyJNu3b8+WLVvyox/9KNdff31uvfXWkdb2JGEYAIBlc6rbLh85ciRXXHFF7rrrrjz66KOZmprK1VdfnVe84hUjP75lEgAAnLHOPffcbNu2LYcPHx7L/MIwAADL5rHHHsvU1FSmpqZy7bXXPu35Y8eO5b777svLXvaysRzfMgkAAJbNqZZJJMkXvvCFbN26Nc973vNy8803C8MAAPTjyTXD42aZBAAA3XJmGACAJMNdCu25RhgGAGDZzM3NPW3fzMxMZmZmJnJ8yyQAAOiWMAwAQLeEYQAAuiUMAwDQLWEYAIBuCcMAAHTLpdUAAJj38RrtfH9r4esWr1q1KpdddtlT23feeWeOHDmS2dnZ7N27N1dffXWS5A1veENuuummkV9yTRgGAGDZnHPOOTl48OBJ+44cOZJNmzblfe9731NheFwskwAA4Ixz+eWXZ/369bnnnnvGehxhGACAZfPYY49lamoqU1NTufbaa0967pZbbsl73/vesR7fMgkAAJbNqZZJPOk1r3lNkuSLX/zi2I7vzDAAAGescZ8dFoYBADhjvfa1r81DDz2U+++/fyzzWyYBAMC8IS6FthxuueWWbN++fSxzDx2Gq2pVkv1JvtNae0NVbUmyJ8mGJAeS/FJr7c/HUiUAAM9Jc3NzT9s3MzNz0vWEr7nmmrQ2nqC+mGUS70xy6ITtDya5rbV2SZKHkuwYZWEAADBuQ4XhqtqU5G8k+e3BdiW5MsknB0N2J3njOAoEAIBxGfbM8K8n+YdJ/mKwvSHJD1prxwfbR5NcPOLaAABgrGqh9RdV9YYkr2+t/b2qmklyU5K3JrlvsEQiVfXiJHe31l5+itfvTLIzSTZu3Lhtz549I/0DnGhubi5r164d2/wjd+DA4sZv2zaeOlaoFddvlkzP+6LffdHvky02IiSnFxPWr1+fSy65ZPEvXKInnngiq1atGslchw8fzsMPP3zSvtnZ2QOttelhXj/MB+heneSaqnp9kjVJzkvyG0nOr6rVg7PDm5J851Qvbq3dnuT2JJmenm4nLoYetX379mWc84/c7Ozixo9p4fhKteL6zZLpeV/0uy/6fbLFRoTk9GLCoUOHsm7dusW/cIkeeeSRkR13zZo12bp162m/fsFlEq21d7fWNrXWNie5Psnvt9b+dpLPJXnTYNgNST5z2lUAAMAyWMpNN341ya9U1eHMryG+YzQlAQCwLKpG+7WAY8eOZWpqKlNTU3nRi16Uiy+++KntqsqNN9741NgPfehD2bVr18j/yIu66UZrbV+SfYPH30ryypFXBABAFzZs2JCDBw8mSXbt2pW1a9fmpptuSjK//OFTn/pU3v3ud+fCCy8cWw1uxwwAwBln9erV2blzZ2677baxHkcYBgDgjPS2t70tH/vYx552tYhREoYBADgjnXfeeXnLW96SD3/4w2M7hjAMAMAZ613velfuuOOOPProo2OZXxgGAOCMdcEFF+S6667LHXeM58JlwjAAAPNaG+3XiNx44435/ve/P7L5TrSoS6sBAMA4/OQ1hOfm5p56vHHjxvzwhz8cy3GdGQYAoFvCMAAA3RKGAQA61ka4tnfSRlG7MAwA0Kk1a9bk2LFjKzIQt9Zy7NixrFmzZknz+AAdAECnNm3alKNHj+Z73/veRI/7+OOPLznEJvNhftOmTUuaQxgGAOjUWWedlS1btkz8uPv27cvWrVsnftxTsUwCAIBuCcMAAHRLGAYAoFvCMAAA3RKGAQDoljAMAEC3hGEAALolDAMA0C1hGACAbgnDAAB0SxgGAKBbwjAAAN0ShgEA6JYwDABAt4RhAAC6JQwDANAtYRgAgG4JwwAAdEsYBgCgW8IwAADdWjAMV9WaqvrvVfWVqvpaVb1nsH9LVX2pqg5X1e9W1U+Nv1wAABidYc4M/yjJla21y5NMJXldVb0qyQeT3NZauyTJQ0l2jK9MAAAYvQXDcJs3N9g8a/DVklyZ5JOD/buTvHEsFQIAwJgMtWa4qlZV1cEkDya5J8n/SvKD1trxwZCjSS4eT4kAADAe1VobfnDV+Uk+neSfJPnoYIlEqurFSe5urb38FK/ZmWRnkmzcuHHbnj17RlH3Kc3NzWXt2rVjm3/kDhxY3Pht28ZTxwq14vrNkul5X/S7L/p9ssVGhGRlxYRx93t2dvZAa216mLGLCsNJUlX/NMljSX41yYtaa8er6meT7Gqt/fVne+309HTbv3//oo63GPv27cvMzMzY5h+5qsWNX2SvnutWXL9ZMj3vi373Rb9PttiIkKysmDDuflfV0GF4mKtJvHBwRjhVdU6Sn09yKMnnkrxpMOyGJJ85vXIBAGB5rB5izEVJdlfVqsyH50+01u6qqq8n2VNV703y5SR3jLFOAAAYuQXDcGvt/iRbT7H/W0leOY6iAABgEtyBDgCAbgnDAAB0SxgGAKBbwjAAAN0ShgEA6JYwDABAt4RhAAC6JQwDANAtYRgAgG4JwwAAdEsYBgCgW8IwAADdEoYBAOiWMAwAQLeEYQAAuiUMAwDQLWEYAIBuCcMAAHRLGAYAoFvCMAAA3RKGAQDoljAMAEC3hGEAALolDAMA0C1hGACAbgnDAAB0SxgGAKBbwjAAAN0ShgEA6JYwDABAt4RhAAC6JQwDANCtBcNwVb24qj5XVV+vqq9V1TsH+y+oqnuq6puD7y8Yf7kAADA6w5wZPp7kxtbaS5O8KsnbquqlSW5Ocm9r7dIk9w62AQBgxVgwDLfWHmit/Y/B40eSHEpycZLtSXYPhu1O8sZxFQkAAOOwqDXDVbU5ydYkX0qysbX2wOCp7ybZONLKAABgzKq1NtzAqrVJ/muS97XWPlVVP2itnX/C8w+11p62briqdibZmSQbN27ctmfPntFUfgpzc3NZu3bt2OYfuQMHFjd+27bx1LFCrbh+s2R63hf97ot+n2yxESFZWTFh3P2enZ090FqbHmbsUGG4qs5KcleSz7bWfm2w7xtJZlprD1TVRUn2tdZ+5tnmmZ6ebvv37x+mrtOyb9++zMzMjG3+kata3Pgh/+PSixXXb5ZMz/ui333R75MtNiIkKysmjLvfVTV0GB7mahKV5I4kh54MwgN7k9wweHxDks8stlAAAFhOq4cY8+okv5Tkj6rq4GDfP0rygSSfqKodSb6d5LrxlAgAAOOxYBhurX0xyTOdrL9qtOUAAMDkuAMdAADdEoYBAOiWMAwAQLeEYQAAuiUMAwDQLWEYAIBuCcMAAHRLGAYAoFvCMAAA3RKGAQDoljAMAEC3hGEAALolDAMA0C1hGACAbgnDAAB0SxgGAKBbwjAAAN0ShgEA6JYwDABAt4RhAAC6JQwDANAtYRgAgG4JwwAAdEsYBgCgW8IwAADdEoYBAOiWMAwAQLeEYQAAuiUMAwDQLWEYAIBuCcMAAHRLGAYAoFvCMAAA3VowDFfVR6rqwar66gn7Lqiqe6rqm4PvLxhvmQAAMHrDnBn+aJLX/cS+m5Pc21q7NMm9g20AAFhRFgzDrbXPJ/nTn9i9PcnuwePdSd444roAAGDsqrW28KCqzUnuaq29fLD9g9ba+YPHleShJ7dP8dqdSXYmycaNG7ft2bNnNJWfwtzcXNauXTu2+UfuwIHFjd+2bTx1rFArrt8smZ73Rb/7ot8nW2xESFZWTBh3v2dnZw+01qaHGbvkMDzYfqi1tuC64enp6bZ///5h6jot+/bty8zMzNjmH7mqxY0folc9WXH9Zsn0vC/63Rf9PtliI0KysmLCuPtdVUOH4dO9msSfVNVFg4NdlOTB05wHAACWzemG4b1Jbhg8viHJZ0ZTDgAATM4wl1b7nSR/kORnqupoVe1I8oEkP19V30zyc4NtAABYUVYvNKC19uZneOqqEdcCAAAT5Q50AAB0SxgGAKBbwjAAAN0ShgEA6JYwDABAt4RhAAC6JQwDANAtYRgAgG4JwwAAdEsYBgCgW8IwAADdEoYBAOiWMAwAQLeEYQAAuiUMAwDQLWEYAIBuCcMAAHRLGAYAoFvCMAAA3RKGAQDoljAMAEC3hGEAALolDAMA0C1hGACAbgnDAAB0SxgGAKBbwjAAAN0ShgEA6JYwDABAt4RhAAC6JQwDANAtYRgAgG6tXsqLq+p1SX4jyaokv91a+8BIqgIAYFHqPbWo8e3WNqZKVpbTPjNcVauS/GaSX0jy0iRvrqqXjqowAAAYt6Usk3hlksOttW+11v48yZ4k20dTFgAAjN9SwvDFSf7PCdtHB/sAAGBFWNKa4WFU1c4kOwebc1X1jTEe7sIk3x/j/MurFrcWqAPP7X5zKnreF/3ui34v1a7FDa9dy5orxt3vvzLswKWE4e8kefEJ25sG+07SWrs9ye1LOM7Qqmp/a216Esdi+el3f/S8L/rdF/3uy5nU76Usk/jDJJdW1Zaq+qkk1yfZO5qyAABg/E77zHBr7XhVvT3JZzN/abWPtNa+NrLKAABgzJa0Zri19ntJfm9EtYzCRJZjcMbQ7/7oeV/0uy/63Zczpt/VmgsuAwDQJ7djBgCgWysyDFfV66rqG1V1uKpuPsXzZ1fV7w6e/1JVbZ58lYzKEP3+lar6elXdX1X3VtXQl1PhzLNQv08Y9zerqlXVGfFpZE7PMP2uqusGP+Nfq6qPT7pGRmuI3+l/uao+V1VfHvxef/1y1MnSVdVHqurBqvrqMzxfVfXhwd+F+6vqFZOuMVmBYXjI20DvSPJQa+2SJLcl+eBkq2RUhuz3l5NMt9b+apJPJvkXk62SURn2Nu9VtS7JO5N8abIVMkrD9LuqLk3y7iSvbq29LMm7Jl4oIzPkz/g/TvKJ1trWzF+p6l9PtkpG6KNJXvcsz/9CkksHXzuT/JsJ1PQ0Ky4MZ7jbQG9Psnvw+JNJrqpyx4oVasF+t9Y+11r74WDzvsxf85qVadjbvP/zzP8n9/FJFsfIDdPvv5vkN1trDyVJa+3BCdfIaA3T85bkvMHj9Un+7wTrY4Raa59P8qfPMmR7kn/f5t2X5Pyqumgy1f1/KzEMD3Mb6KfGtNaOJ3k4yYaJVMeoLfa23zuS3D3WihinBfs9eBvtxa21/zTJwhiLYX6+fzrJT1fVf6uq+6rq2c4yceYbpue7kvxiVR3N/BWrfnkypbEMFvtv/FiM/XbMMClV9YtJppP8teWuhfGoqucl+bUkb13mUpic1Zl/C3Um8+/6fL6qLmut/WBZq2Kc3pzko621f1VVP5vkP1TVy1trf7HchfHctBLPDA9zG+inxlTV6sy/zXJsItUxakPd9ruqfi7JLUmuaa39aEK1MXoL9Xtdkpcn2VdVR5K8KsleH6JbsYb5+T6aZG9r7cettf+d5H9mPhyzMg3T8x1JPpEkrbU/SLImyYUTqY5JG+rf+HFbiWF4mNtA701yw+Dxm5L8fnNB5ZVqwX5X1dYkv5X5IGw94cr2rP1urT3cWruwtba5tbY582vEr2mt7V+eclmiYX6f35n5s8Kpqgszv2ziW5MskpEapud/nOSqJKmql2Q+DH9volUyKXuTvGVwVYlXJXm4tfbApItYccsknuk20FX1z5Lsb63tTXJH5t9WOZz5hdvXL1/FLMWQ/f6XSdYm+Y+Dz0n+cWvtmmUrmtM2ZL95jhiy359N8tqq+nqSJ5L8g9aad/pWqCF7fmOSf1dVfz/zH6Z7qxNaK1NV/U7m/zN74WAN+K1JzkqS1tq/zfya8NcnOZzkh0n+zrLU6e8XAAC9WonLJAAAYCSEYQAAuiUMAwDQLWEYAIBuCcMAAHRLGAYAoFvCMAAA3RKGAQDo1v8Dq/8WDFGP9fgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "TP, FP, FN, TN = confusion_table(mlp(X, 1.57166185e+01,  1.58250335e+01, -3.52251392e+00, -7.49692230e+00,\n",
    "       -5.89003641e+01,  3.80089003e-02,  2.37483939e-01,  2.87007844e-01,\n",
    "       -9.48869044e-01,  1.39257951e+00, -1.28133971e-01, -5.79757478e-01,\n",
    "        3.16485573e-01, -1.03528732e+00,  1.34439583e+00,  1.12009673e+02,\n",
    "       -5.14928444e+01, -5.51022639e+01, -3.22706767e+01), Y)\n",
    "show_result(TP, FP, FN, TN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "頑張ったわりには、ロジスティック回帰と同じ性能だった(´・ω・`)\n",
    "\n",
    "# 結論\n",
    "\n",
    "多層パーセプトロン (Multilayer perceptron, MLP)を scipy.optimize.curve_fit で実装する試みは、出力層が１つだけなら、可能でした。出力層が２つ以上の場合は、どうすれば良いかわかりません。無理かも。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}