{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"scipy.optimize.curve_fit に関する記事を次々書いているところですが、\n",
"\n",
"* 少ない観測値を補間してから、正規分布の線形和で近似する\n",
"* カーブフィッティング手法 scipy.optimize.curve_fit の使い方を理解する\n",
"* ロジスティック回帰を scipy.optimize.curve_fit で実装する \n",
"\n",
"に続いて、今回は「多層パーセプトロン (Multilayer perceptron, MLP)を scipy.optimize.curve_fit で実装する」に挑戦します。\n",
"\n",
"# 多層パーセプトロン (Multilayer perceptron, MLP) とは\n",
"\n",
"過去に \n",
"\n",
"* 多層パーセプトロン (Multilayer perceptron, MLP)をExcelで理解する\n",
"* 多層パーセプトロン (Multilayer perceptron, MLP)をPythonで理解する\n",
"\n",
"を書いたので、そちらをご覧ください。前者の方が内容は平易だと思います。\n",
"\n",
"# 「あやめのデータ」を読み込む\n",
"\n",
"「あやめのデータ」については 実習用データ をご参照ください。"
]
},
{
"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",
"多層パーセプトロン (Multilayer perceptron, MLP)をExcelで理解する と同じく、入力層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\u001b[0m in \u001b[0;36m\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": [
"データも取り直し。ロジスティック回帰を scipy.optimize.curve_fit で実装する と同じものです。"
]
},
{
"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": [
"結果表示の関数も ロジスティック回帰を scipy.optimize.curve_fit で実装する と同じ。"
]
},
{
"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": [
""
]
},
"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
}