{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"オリジナル作成: 2011/06/04"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\t
\n", "\t\t第1章-Sageを使って線形回帰を試してみるでは、\n", "\t\t与えられたモデルパラメータM, α、βに対して解いていましたが、\n", "\t\t実際には自分でベストなモデルパラメータM, α、βを求める必要があります。\n", "\t\t\n", "\t\tパターン認識と機械学習\n", "\t\tの3章ではエビデンス関数を使って最適なパラメータを求める方法が説明されています。\n", "\t
\t\n", "\t\n", "\t\tここでは、Sageを使ってエビデンス関数の評価、最適なモデルパラメータの推定を試してみます。\n", "\t
\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t\n", "\t\tエビデンス関数の対数表現は、式(3.86)\n", "$$\n", "\t\t\\ln p(t|\\alpha, \\beta) = \\frac{M}{2} \\ln \\alpha + \\frac{N}{2} \\ln \\beta -E(m_N) - \\frac{1}{2} \\ln | A | + \\frac{N}{2} \\ln (2 \\pi)\n", "$$\n", "\t\tで与えられます。\n", "\t
\n", "\t\n", "\t\tここで、$E(m_N)$は、式(3.82)\n", "$$\n", "\t\tE(m_N) = \\frac{\\beta}{2}|| t = \\Phi m_N ||^2 + \\frac{\\alpha}{2} m_N^T m_N\n", "$$\n", "\t\t$m_N$は、式(3.84)\n", "$$\n", "\t\tm_N = \\beta A^{-1} \\Phi^T t\n", "$$\n", "\t\tAは、式(3.81)\n", "$$\n", "\t\tA = \\alpha I + \\beta \\Phi^T \\Phi\n", "$$\n", "\t
\n", "\t\n", "\t\t目指すは、図3.14です。\t\t\n", "\t
\n", "\t\n", "\t\t\n", "\t
\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t\n", "\t\t第1章-Sageを使って線形回帰を試してみる\n", "\t\tで使ったデータと同じものを座標Xと目的値tにセットし、関数Φを定義します。\n", "\t
\n", "" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# PRML fig.3.14の再現\n", "# PRMLのsin曲線のデータ\n", "data = matrix([\n", " [0.000000, 0.349486],\n", " [0.111111, 0.830839],\n", " [0.222222, 1.007332],\n", " [0.333333, 0.971507],\n", " [0.444444, 0.133066],\n", " [0.555556, 0.166823],\n", " [0.666667, -0.848307],\n", " [0.777778, -0.445686],\n", " [0.888889, -0.563567],\n", " [1.000000, 0.261502],\n", " ]);\n", "X = data.column(0)\n", "t = data.column(1)\n", "# データを増やす場合\n", "# N = 25\n", "# X = vector([random() for i in range(25)])\n", "# t = vector([(sin(2*pi*x) + +gauss(0, 0.2)).n() for x in X.list()])" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Φ関数定義\n", "def _phi(x, j):\n", " return x^j" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t\n", "\t\tまずは、α、βを固定値\n", "\t\t
\n", "\t\t一番良いエビデンスは、M=4あたりになります。\n", "\t\t理由は分かりませんが、図3.14のようなM=3以降の急激な下降は見られません。\n", "\t
\n", "" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEdCAYAAADkeGc2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAF0ZJREFUeJzt3XtQFef9x/HPclEDHBQigmKi0MSqRK2aQRoco9Gpl0Sp\nUROO0do0Gls0Rs1UbTqtvcQxOnZspyaOYzpNOxMkaa0TL3FMg6Rao9QLYCZRbKvBekFzASHoDwX2\n9wcRg6AI5xx24Xm/Zphh9+x5nq87zuc8PLv7HMu2bVsAAKMEOV0AAKD1Ef4AYCDCHwAMRPgDgIEI\nfwAwEOEPAAYi/AHAQIQ/ABiI8AcAAxH+AGAgV4b/K6+8ooSEBN11111KSUnRwYMHnS7JNVauXKnk\n5GRFRkYqNjZWkydP1okTJ5wuy7VWrlypoKAgLV682OlSXOfcuXOaOXOmunbtqrCwMA0aNEhHjhxx\nuizXqKmp0c9+9jMlJiYqLCxM9913n1566SWny/Ib14X/m2++qRdeeEG//OUvlZeXp0GDBmns2LH6\n7LPPnC7NFfbu3avnnntOubm5eu+993Tt2jV95zvf0ZUrV5wuzXUOHjyojRs3atCgQU6X4jqlpaVK\nTU1Vx44dtWvXLh07dky/+c1vFBUV5XRprvHyyy9rw4YNevXVV3X8+HGtXr1aq1ev1rp165wuzT9s\nlxk2bJi9YMGCuu2amho7Pj7eXrVqlYNVudenn35qW5Zl79271+lSXKW8vNzu06ePnZ2dbY8cOdJe\ntGiR0yW5ytKlS+0RI0Y4XYarPfbYY/bs2bPr7ZsyZYo9c+ZMhyryL1eN/K9du6bDhw9r9OjRdfss\ny9KYMWO0f/9+Bytzr9LSUlmWpejoaKdLcZV58+Zp4sSJeuSRR5wuxZW2bdumBx98UE888YRiY2M1\nZMgQvfbaa06X5SoPPfSQsrOz9e9//1uSVFBQoH379mnChAkOV+YfIU4X8HWfffaZqqurFRsbW29/\nbGysCgsLHarKvWzb1sKFCzV8+HD179/f6XJcIysrS/n5+Tp06JDTpbjWyZMntX79er3wwgv66U9/\nqtzcXC1YsECdOnXSjBkznC7PFZYtW6aysjL17dtXwcHBqqmp0YoVK5Senu50aX7hqvC/Fdu2ZVmW\n02W4TkZGhj7++GPt27fP6VJc48yZM1q4cKH+/ve/KzQ01OlyXKumpkbJycn69a9/LUkaNGiQPvro\nI61fv57w/8qbb76pzMxMZWVlqX///srPz9fzzz+vHj16aObMmU6X5zNXhX/Xrl0VHBysCxcu1Nt/\n8eLFBn8NmG7+/Pl65513tHfvXnXv3t3pclzj8OHD+vTTTzV06FDZX31PUXV1tfbs2aN169apsrKS\ngYSk7t27q1+/fvX29evXT3/7298cqsh9lixZohdffFHTpk2TJCUlJemTTz7RypUr20X4u2rOPzQ0\nVEOHDlV2dnbdPtu2lZ2drYceesjBytxl/vz5evvtt5WTk6N7773X6XJcZcyYMfrwww+Vn5+vgoIC\nFRQU6MEHH9SMGTNUUFBA8H8lNTW1wVRqYWGhevXq5VBF7nP58uUG/1+CgoJUU1PjUEX+5aqRvyQt\nXrxYs2bN0tChQ5WcnKy1a9fq8uXL+v73v+90aa6QkZGhTZs2aevWrQoPD6/7K6lz587q1KmTw9U5\nLzw8vMH1j/DwcN19990NRromW7RokVJTU7Vy5Uo98cQTys3N1WuvvaaNGzc6XZprTJw4UStWrNA9\n99yjpKQkHTlyRGvXrtXs2bOdLs0/nL3ZqHGvvPKK3atXL7tTp052SkqKffDgQadLcg3LsuygoKAG\nP3/605+cLs21Ro0axa2ejdixY4c9YMAA+6677rL79+9v/+EPf3C6JFf58ssv7UWLFtm9e/e2w8LC\n7Pvuu8/++c9/bl+7ds3p0vzCsm2+wB0ATOOqOX8AQOsg/AHAQIQ/ABiI8AcAAxH+AGAgwh8ADOTY\nQ162bau8vNyp7gGg3fJ4PE0+ze5Y+JeXl6tz585OdQ8A7dalS5cUGRl522Mce8irqZF/enq6srKy\nfO7HX+34s6322o4/22qv7fizrfbajj/baq/tNNWWq0f+lmXd9pMpJCSkyU+uO+GvdvzZVnttx59t\ntdd2/NlWe23Hn22113b80Va7v+Dr9Xpd15Y/a/IHzlHTOEdN4xw1zU3nyLVr+0yaNElbt251ugxX\n4xw1jXPUNM5R09rjOWr3I38AQEPBv/jFL37hdBG3MmDAAKdLcD3OUdM4R03jHDWtvZ0j1077AAAC\nh2kfADAQ4Q8ABiL8AcBAhD8AGIjwBwADtUr4P/300woKCqr3M2HChFseb9u2ysrKxI1IABAYrba2\nz/jx4/X666/XBXrHjh1veez1FT/vZGU6AEDztVr4d+zYUTExMa3VHQDgNlptzv/9999XbGys+vbt\nq4yMDH3xxRet1TUA4Cat8oTvW2+9pbCwMCUkJOi///2vfvKTn8jj8Wj//v2NrjldVlbGtA8ABJDf\nwz8zM1Nz586tbdyytHPnTqWmptY75tSpU/rGN76h7OxsjRo1qkEb18N//PjxCgmpPzPl9Xpdt0wr\nALQ1fg//iooKXbhwoW47Pj6+0Yu73bp104oVKzRnzpwGrzHyB4DA8vsF3/DwcCUmJt72mDNnzujz\nzz9X9+7d/d09AOAOBPyCb0VFhZYsWaLc3FwVFRUpOztb3/3ud9WnTx+NHTs20N0DABoR8Fs9g4OD\ndfToUf35z39WaWmpevToobFjx+pXv/qVQkNDA909AKARrlzPnzl/AAgs1vYBAAMR/gBgIMIfAAxE\n+AOAgQh/ADCQq8M/PT1dkyZN0qZNm5wuBQDaFW71BAADuXrkDwAIDMIfAAxE+AOAgQh/ADAQ4Q8A\nBiL8AcBAhD8AGIjwBwADEf4AYCBXhz/LOwBAYLC8AwAYyNUjfwBAYBD+AGAgwh8ADET4A4CBCH8A\nMBDhDwAGIvwBwECEPwAYyNXhzxO+ABAYPOELAAZy9cgfABAYhD8AGIjwBwADEf4AYCDCHwAMRPgD\ngIEIfwAwEOEPAAYi/AHAQK4Of5Z3AIDAYHkHADCQq0f+AIDAIPwBwECEPwAYyOfw37Jli8aNG6eY\nmBgFBQXp6NGjDY6prKzUvHnz1LVrV3k8Hk2dOlUXL170tWsAQAv5HP4VFRUaPny4Vq1aJcuyGj1m\n4cKF2rFjhzZv3qw9e/bo3LlzmjJliq9dAwBayG93+xQVFSkhIUH5+fkaOHBg3f6ysjLFxMQoKytL\nkydPliQVFhaqX79+OnDggJKTkxu0xd0+ABBYAZ/zP3z4sKqqqjR69Oi6fd/85jd17733av/+/YHu\nHgDQiICHf3FxsTp06NBgBB8bG6vi4uJAdw8AaESzwj8zM1Mej0cej0eRkZHat29fizu2bfuW1wgA\nAIEV0pyD09LSlJKSUrcdHx/f5Hvi4uJ09epVlZWV1Rv9X7x4UbGxsbd9b3p6ukJC6pfo9Xrl9Xqb\nUzYA4CbNCv/w8HAlJibe8vXGRvJDhw5VSEiIsrOz6y74njhxQqdPn9a3v/3t2/aXlZXFBV8ACIBm\nhX9jSkpKdPr0aZ09e1a2bev48eOybVtxcXGKjY1VZGSknnnmGS1evFhRUVHyeDxasGCBUlNTG73T\nBwAQeD5f8N26dasGDx6siRMnyrIseb1eDRkyRBs2bKg7Zu3atXrsscc0depUjRw5Uj169NDmzZt9\n7RoA0EKs6gkABmJtHwAwEOEPAAYi/AHAQIQ/2q3du6Vt26TKSqcrAdyH8Ee79KMfSaNHS5MmSY88\nIl296nRFgLu4+m6f8ePHKyQkhKd60SwVFVJERP19778vPfywI+UAruTzQ16BxBO+bcOVK9KTT0q7\ndkkDB0pbtkg9ezpXT4cOUnh47YfAdVFRztUDuBHTPvDZ2rW1c+tXr0qHDkkvvOBsPaGh0htv1AZ+\nx47SihW1H0oAbnD1yB9tw6ef1t/+/HNn6vi6tDTpiy8k25ZYPBZoiJE/fPaDH0jXZ+eCg6WMDGfr\n+TqCH2gcI3/4bMAA6ehRaf9+qX9/plhuJS9Pevdd6YEHpEcfdboamI7wh1/06lX7g8YdOFB7t9H1\nW05//3tp/nxna4LZmPYBWsHmzfWfNcjKcq4WQCL8gVZx83cgJSQ4UwdwHdM+QCuYO1cqLJR27Kid\n8//tb52uCKZz9RO+rOcPAIHh6mmf9PR0TZo0SZs2bXK6FABoVxj5A4Y6ckSaN08qL5eWLpVmznS6\nIrQmwh8wVHy8dO5c7e/BwdKHH0r9+jlbE1qPq6d9AATG//3fjeCXpOpqqajIuXrQ+gh/wECdOtWu\nf3Rdr15SSopz9aD1casnYKi33pJef10qK5NmzJC6dHG6IrQm5vwBuMaePbXPQvTtKz39tNPVtG+M\n/AG4wj//WfuVm9XVtdunT0vLlztbU3vGnD8AV3jnnRvBL9V+QRACh/AH4ApJSbffhn+5etonPT2d\nL3AHDPHUU7VTPW+/Xfu8gVvWP7Jtqaqq9utB2xMu+ALALWzfXnsnVHm5tGiRtGaN0xX5D+EPALcQ\nFSWVlt7Y3rtXGj7cuXr8iTl/AGhEdbX05Zf191265EwtgUD4A0AjgoOlJUtubCcnS6NHO1ePvxH+\nAHALK1ZI+/dLO3dK//hH7bIYTqqulp55pvZp7ORk6dSplrfFnD8AtBEbN0rPPntje9y42g+mlmDk\nDwBtRHHx7bebg/AHgDbiySfrL8A3d27L23L1Q14AgBv69JHy86Xdu6X77/fttlPm/AHAQK6e9uEL\n3AEgMBj5A4CBXD3yBwAEBuEPAAYi/AHAQD6H/5YtWzRu3DjFxMQoKChIR48ebXDMyJEjFRQUVPcT\nHBysjIwMX7sGALSQz+FfUVGh4cOHa9WqVbIsq9FjLMvSs88+qwsXLqi4uFjnz5/X6tWrfe0aANBC\nPj/kNWPGDElSUVGRbnfjUFhYmGJiYnztDgDgB6025//GG28oJiZGAwYM0IsvvqgrV660VtcAgJu0\nyvIOTz31lHr16qUePXro6NGjWrJkiU6cOKG//vWvrdE9AOAmzQr/zMxMzf1qJSHLsrRz506lpqY2\n+b7Zs2fX/Z6UlKS4uDiNGTNGp06dUkJCQjNLdkZNjRTEvVEA2olmhX9aWppSUlLqtuPj41vU6bBh\nw2Tbtv7zn//cNvzT09MVElK/RK/XK6/X26J+W6K8XJo8uXYhpUGDpG3bpJ49W617AAiIZoV/eHi4\nEhMTb/n6re72uVleXp4sy1L37t1ve1xWVpbjyzusXi1lZ9f+np8vLV0qvfGGoyUBgM98nvMvKSnR\n6dOndfbsWdm2rePHj8u2bcXFxSk2NlYnT55UZmamJkyYoLvvvlsFBQVavHixHn74YT3wwAP++DcE\nVEnJ7bcBoC3yeRZ769atGjx4sCZOnCjLsuT1ejVkyBBt2LBBktShQwe99957Gjt2rPr166cf//jH\nmjZtmrZu3epz8a3h2Welzp1rf+/QQVqwwNl6AMAfWNXzDpw5Ix08KCUl1X6ZAgC0dXyT1x3o2ZOL\nvADaF25eBAADEf4AYCDCHwAMRPgDgIFcHf58gTsABAa3egKAgVw98gcABAbhDwAGIvwBwECEPwAY\niPAHAAMR/gBgIMIfAAxE+AOAgQh/ADCQq8Of5R0AIDBY3gEADOTqkT8AIDAIfwAwEOHfRhUUSLm5\nUk2N05UAaIsI/zZo6VLpW9+SUlKkxx/nAwBA83HBt425dEnq0qX+vgMHpGHDnKkHQNvEyL+NCQ2V\nQkLq77vrLmdqAdB2Ef5tTFiYtH79jQ+AJUukgQOdrQlA28O0Txt1+bJUVSVxegC0REjTh8CNwsKc\nrgBAW+bqaR+WdwCAwGDaBwAM5OqRPwAgMAh/ADAQ4Q8ABiL8AcBAhD8AGIjwBwADEf4AYCDCHwAM\n5Orw5wlfAAgMnvAFAAO5euQPAAgMwh8ADET4A4CBfAr/qqoqLV26VAMHDlRERITi4+M1a9YsnT9/\nvt5xJSUleuqpp9S5c2dFRUVp9uzZqqio8KlwAEDL+RT+ly9fVn5+vpYvX668vDxt2bJFhYWFSktL\nq3fc9OnTdezYMWVnZ2vHjh3as2eP5s6d61PhAICW8/vdPocOHdKwYcNUVFSknj176tixY0pKStLh\nw4c1ePBgSdKuXbv06KOP6syZM4qLi2vQBnf7AEBg+X3Ov7S0VJZlqUuXLpKkAwcOKCoqqi74JWnM\nmDGyLEu5ubn+7h4AcAf8Gv6VlZVatmyZpk+froiICElScXGxunXrVu+44OBgRUdHq7i42J/dAwDu\nULPCPzMzUx6PRx6PR5GRkdq3b1/da1VVVZo2bZosy9Krr77aZFu2bcuyrOZXDADwWUhzDk5LS1NK\nSkrddnx8vKQbwf+///1Pu3fvrhv1S1JcXJwuXrxYr53q6mqVlJQoNjb2tv2lp6crJKR+iV6vV16v\ntzllAwBu4vMF3+vBf/LkSeXk5Cg6Orre68ePH1dSUpIOHTpUN+//7rvvasKECVzwBQCHNGvkf7Pq\n6mpNmTJF+fn52r59u65du6YLFy5IkqKjoxUaGqq+fftq7NixmjNnjtavX6+rV6/queeek9frbTT4\nv/hCWr269vezZyWyHwD8z6eRf1FRkRITE+vtuz6Xn5OToxEjRkiqvQNo/vz52rZtm4KCgjR16lT9\n7ne/U1hYWL331tRIDz4o5eWVSeqse+65pI8+ipTH09IKAQCNcdWqnmfPSj17SlJt+EuXlJsbqeRk\nZ+sCgPbGVWv7xMRIX58JCg+Xevd2rBwAaLdcFf4dOki7dkljxtRu/+Uv0k2PCAAA/MBV0z7XcbcP\nAASWq0b+AIDWQfgDgIFcHf58gTsABAZz/gBgIFeP/AEAgUH4A4CBCH8AMBDhDwAGIvwBwECEPwAY\niPAHAAMR/gBgIMIfAAzk6vBneQcACAyWdwAAA7l65A8ACAzCHwAMRPgDgIEIfwAwEOEPAAYi/AHA\nQIQ/ABiI8AcAAxH+AGAgV4c/yzsAQGCwvAMAGMjVI38AQGAQ/gBgIMIfAAxE+AOAgQh/ADAQ4Q8A\nBiL8AcBAhD8AGMjV4c8TvgAQGDzhCwAGcvXIHwAQGIQ/ABiI8AcAA/kU/lVVVVq6dKkGDhyoiIgI\nxcfHa9asWTp//ny943r37q2goKC6n+DgYK1evdqnwgEALRfiy5svX76s/Px8LV++XAMHDlRJSYkW\nLFigtLQ0/etf/6o7zrIsvfTSS5ozZ46uX1/2eDy+VQ4AaDGfwj8yMlK7du2qt2/dunUaNmyYzpw5\no549e9btj4iIUExMjC/dAQD8xO9z/qWlpbIsS126dKm3/+WXX1bXrl01ZMgQrVmzRtXV1f7uGgBw\nh3wa+d+ssrJSy5Yt0/Tp0xUREVG3//nnn9eQIUMUHR2tDz74QMuWLVNxcbHWrFnjz+4BAHeoWQ95\nZWZmau7cubVvtCzt3LlTqampkmov/j7++OM6f/68cnJy6oX/zf74xz/qhz/8ob788kuFhoY2eJ2H\nvAAgsJoV/hUVFbpw4ULddnx8vDp27KiqqipNmzZNn3zyiXbv3q2oqKjbtvPxxx9rwIABOn78uO6/\n//4Gr18P//HjxyskpP4fJ16vV16v905LBgA0olnTPuHh4UpMTKy373rwnzx5Ujk5OU0GvyTl5eUp\nKChI3bp1u+1xWVlZjPwBIAB8mvOvrq7WlClTlJ+fr+3bt+vatWt1fxlER0crNDRUBw4cUG5urkaN\nGiWPx6MPPvhAixcv1syZM9W5c2e//CMAAM3j08JuRUVFDf4SsG1blmUpJydHI0aMUF5enjIyMlRY\nWKjKykolJCToe9/7nhYtWtTofL/EnD8ABBqregKAgVjbBwAMRPgDgIEIfwAwEOEPAAYi/AHAQIQ/\nABjI1eGfnp6uSZMmadOmTU6XAgDtCvf5A4CBXD3yBwAEBuEPAAYi/AHAQK6c87dtW+Xl5fJ4PLIs\ny+lyAKDdcWX4AwACi2kfADAQ4Q8ABiL8AcBAhD8AGIjwBwADEf4AYCDCHwAM9P+5xtKo/VO/4QAA\nAABJRU5ErkJggg==\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# α=5*10^-3, β=11.1を固定して計算すると\n", "ln_p_List = []\n", "N = len(t)\n", "alpha = 5*10^-3\n", "beta = 11.1\n", "for M in range(10):\n", " Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1))] for x in X.list()])\n", " Phi_t = Phi.transpose()\n", " # m_N, Aを定義\n", " A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi\n", " m_N = beta*A.inverse() * Phi_t * t\n", " # エビデンスの対数\n", " res = (t - Phi*m_N)\n", " E_mN = beta/2* res * res + alpha/2 * m_N * m_N\n", " ln_p_t = M/2 * ln(alpha) + N/2 * ln(beta) - E_mN - 1/2 * ln(A.det()) - (N/2 * ln(2*pi))\n", " ln_p_List += [ln_p_t]\n", "list_plot(ln_p_List, ymin = -26, ymax = -5).show(figsize=4)\n", "# 結果は、M=3以降でfig.3.14のような急激な降下は見られない" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t\n", "\t\tそこで、βを平均値の分散$\\beta_{ML}$を使って計算してみました。\n", "$$\n", "\t\t\\frac{1}{\\beta_{ML}} = \\frac{1}{N} \\sum (y(x_n, W_ML) - t_n)^2 \n", "$$\n", "\t
\n", "\t\n", "\t\t期待に反し、M=0でも高い値となり、あまり良くありません。\n", "\t
\n", "" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEcCAYAAAAvJLSTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAFl1JREFUeJzt3XtwFeX9x/HPniQCuZKQECQQSMRKoEChFKihKhXl0mIU\nwUm4FKxYWqBA6FSwrTq1MCDS6WU6pg60RWxDenGwIGVQY5AOQmyQQOUqxQYpJLSYC4QaDNnfH+dH\n6CGQEM7Z7JLn/Zo549nNnn2+WeInT57dZ9eybdsWAMAoPrcLAAC0PcIfAAxE+AOAgQh/ADAQ4Q8A\nBiL8AcBAhD8AGIjwBwADeTL8bdvWr3/9azH/rHnr1693uwTP4xi1jGPUsvZ4jDwZ/mfPntVjjz2m\ns2fPul2Kp7XHH8hQ4xi1jGPUsvZ4jDwZ/gAAZ7X78A/lb+xQ7ctrvQiOUcs4Ri3jGLXMS8eI8Hdh\nX/xAtt1+QoVj1DKOUcu8dIzCQ1RHq9m2fc0x/ZqamoD/BqO+vj4k+wnlvtrrfkK5r/a6n1Duq73u\nJ5T7aq/7aWlfMTExsiyr2c9bbt3SuaamRnFxcW40DQDtWnV1tWJjY5vdxrXwb6nn37NnT3300Uct\nfgMAgEDX0/N3bdjHsqwWgz02NpbwBwAHtPsTvgCApgh/ADAQ4Q8ABiL8AcBAhP9NZvny5Ro2bJhi\nY2OVnJyshx56SEeOHHG7LM9avny5fD6fFi1a5HYpnnPy5ElNnz5diYmJioyM1KBBg/Tee++5XZZn\nNDQ06KmnnlJ6eroiIyPVp08fLV261O2yQsa1q31wY/7617/q29/+toYOHar6+no9+eSTuv/++3Xw\n4EF16tTJ7fI85W9/+5tWr16tQYMGuV2K51RVVSkzM1P33nuvtm7dqsTERH3wwQeKj493uzTPWLFi\nhV588UWtW7dO/fr1U0lJiWbOnKnOnTtr3rx5bpcXNML/JvOXv/wlYHnt2rXq2rWrdu/erZEjR7pU\nlfecO3dO06ZN05o1a/SjH/3I7XI8Z8WKFUpNTdWaNWsa1/Xq1cvFirxn586dysrK0tixYyVJqamp\nys/P17vvvutyZaHh6WGf7OxsPfDAA567P4eXVFVVybIsJSQkuF2Kp8ydO1cTJkzQl7/8ZbdL8aRN\nmzZp6NCheuSRR5ScnKwhQ4YE/CKAdOedd6qwsFAffPCBJGnv3r3asWOHxo8f73JloeHpnn9BQQGT\nvJph27YWLlyokSNHql+/fm6X4xkFBQUqLS1VSUmJ26V41rFjx5SXl6fvfOc7+v73v6/i4mLNnz9f\nHTt21LRp09wuzxOWLFmimpoa9e3bV2FhYWpoaNCyZcuUnZ3tdmkh4enwR/PmzJmjAwcOaMeOHW6X\n4hknTpzQwoUL9cYbbygiIsLtcjyroaFBw4YNaxwSGzRokPbv36+8vDzC///9/ve/V35+vgoKCtSv\nXz+VlpZqwYIF6t69u6ZPn+52ecGzPai6utqWZFdXV7tdimfNnTvXTk1NtcvKytwuxVNeffVV2+fz\n2REREXZ4eLgdHh5uW5bVuK6hocHtEj2hV69e9uOPPx6wLi8vz+7Ro4dLFXlPz5497by8vIB1S5cu\ntTMyMlyqKLTo+d+E5s2bpz//+c96++23lZqa6nY5njJ69Gj9/e9/D1g3c+ZMZWRkaMmSJS3e7MoU\nmZmZOnz4cMC6w4cPc9L3f5w/f77Jz4vP51NDQ4NLFYUW4X+TmTNnjtavX6+NGzcqKipKFRUVkqS4\nuDh17NjR5ercFxUV1eT8R1RUlLp06aKMjAyXqvKe3NxcZWZmavny5XrkkUdUXFysNWvWaPXq1W6X\n5hkTJkzQsmXL1LNnT/Xv31/vvfeefvKTn2jWrFlulxYabv/pcTUM+1zbpSGMK18vvfSS26V51qhR\no+zc3Fy3y/CczZs32wMGDLA7depk9+vXz/7Vr37ldkmecu7cOTs3N9fu3bu3HRkZaffp08d++umn\n7U8//dTt0kLCtfv5N+fSg16u54EEAIDW8/R1/gAAZxD+AGAgT4c/M3wBwBmM+QOAgTzd8wcAOIPw\nBwADEf4AYCDCHwAMRPgDgIEIfwAwEOEPAAbydPgzyQsAnMEkLwAwkKd7/gAAZxD+AGAgwh8ADET4\nA4CBCH8AMBDhDwAGIvwBwECEPwAYyNPhzwxfAHAGM3wBwECe7vkDAJxB+AOAgdok/B999FH5fL6A\n1/jx49uiaQDAVYS3VUPjxo3T2rVrdekUQ4cOHdqqaQDAFdos/Dt06KCkpKS2ag4A0Iw2G/Pftm2b\nkpOT1bdvX82ZM0cff/xxWzUNALhCm1zq+Yc//EGRkZFKS0vTP/7xDz355JOKiYnRzp07ZVlWk+25\n1BMAnBXy8M/Pz9fs2bP9O7csbdmyRZmZmQHbfPjhh7rttttUWFioUaNGNdkH4Q8Azgr5mH9WVpZG\njBjRuJySktJkm7S0NCUmJuro0aNXDf9LsrOzFR4eWGJOTo5ycnJCVzAAGCjk4R8VFaX09PRmtzlx\n4oTOnDmjW2+9tdntCgoK6PkDgAMcP+FbW1urJ554QsXFxSorK1NhYaEefPBBfeYzn9GYMWOcbh4A\ncBWOX+oZFhamffv2ad26daqqqlL37t01ZswYPfvss4qIiHC6eQDAVXBjNwAwEPf2AQADEf4AYCDC\nHwAMRPgDgIEIfwAwkKfDn2f4AoAzuNQTAAzk6Z4/AMAZhD8AGIjwBwADEf4AYCDCHwAMRPgDgIEI\nfwAwkKfDn0leAOAMJnkBgIE83fMHADiD8AcAAxH+AGAgwh8ADET4A4CBCH8AMBDhDwAGIvwBwECe\nDn9m+AKAM5jhCwAG8nTPHwDgDMIfAAxE+AOAgQh/ADAQ4Q8ABiL8AcBAhD8AGIjwBwADeTr8meEL\nAM5ghi8AGMjTPX8AgDMIfwAwUNDhv2HDBo0dO1ZJSUny+Xzat29fk23q6uo0d+5cJSYmKiYmRpMm\nTdLp06eDbRq4puPHpRkzpEmTpJ073a4G8J6gw7+2tlYjR47Uc889J8uyrrrNwoULtXnzZr3yyiva\nvn27Tp48qYcffjjYpoFruv9+ad066ZVXpDFjpJMn3a4I8JaQnfAtKytTWlqaSktLNXDgwMb1NTU1\nSkpKUkFBgR566CFJ0uHDh5WRkaFdu3Zp2LBhTfbFCV8Eo7pa6tw5cN2bb0r33utOPYAXOT7mv3v3\nbtXX1+ve//k/74477lBqaqp28vc4HBAXJ33uc5eXExKkAQPcqwfwonCnGygvL9ctt9zSpAefnJys\n8vJyp5uHoV5/XVq2TKqtlebPl7p2dbsiwFta1fPPz89XTEyMYmJiFBsbqx07dtxww7ZtX/McAW4+\n+fn+E6w//rF08aLb1UhJSdJPfyqtXk2vH7iaVvX8s7KyNGLEiMbllJSUFj/TrVs3XbhwQTU1NQG9\n/9OnTys5ObnZz2ZnZys8PLDEnJwc5eTktKZsOOyVV6SpUy8vV1dLzz7rXj0AWtaq8I+KilJ6evo1\nv361nvznP/95hYeHq7CwsPGE75EjR3T8+HF98YtfbLa9goICTvjeBN5+O3B5+3Z36gBw/YI+4VtZ\nWam9e/dq//79sm1bhw4d0t69e1VRUSFJio2N1WOPPaZFixZp27Zt2r17tx599FFlZmZe9Uof3HyG\nDw9c5p+1qU8+kSZOlDp0kIYMkcrK3K4IxrODtHbtWtuyLNvn8wW8fvjDHzZu88knn9jz5s2zu3Tp\nYkdHR9uTJk2yKyoqrrnP6upqW5JdXV0dbHloIy+8YNtZWbb91FO2feGC29V4z/PP27Z0+fXgg25X\nBNNxYzegDSxeLK1ceXn5S19ieAzu4t4+QBuYOfPyxLOwMGnuXFfLAZy/zh+AlJEh7dsn7dgh9e0b\nOAkNcAPDPoDBtm+Xzp713/qiY0e3q0FbYtgHMFRurnT33dJXvyrdc4//iiSYw9M9/3Hjxik8PJyJ\nXUCI1dVJnTr5rz265PXXpfvuc68mtC1Pj/kzyQtwRkSEFB3tH/K5JD7evXrQ9hj2AQzk8/nvxxQf\n7/9F8PTT0tChblflf+7CH/8o7d3rdiXtn6eHfTjhCzivocH/y8BtR49KI0ZIZ87463n5ZWnKFLer\nar888E/e1L//7XYFgDm8EPyS/8lrZ8743zc0SD//ubv1tHce+We/7PnnpT59/O8XL3a3FgBtp0uX\nwOWEBHfqMIWnhn0qK/0/ALZdIylOUrXefz9W/fu7XRkAp124IOXkSJs2SXfcIb36qnTbbe7WVFHh\nPx9SVSXNm+e/LUd74anwP3NGSkyUpMvhX1oaq0GD3K0LgJm+8AWppMT/PjJSev99KS3N3ZpCxVPD\nPl26SD/4weXlqVNF8ANwRX395eCXpPPn/bfoaC881fO/pLS0RoMHc7UPAHfdeae0c6f/fXS0tH+/\nlJrqbk2h4qme/yWXHhaWnZ2tBx54QOvXr3e3IABG2rRJWrDA/3zqwkL3g7++Xpo+XYqK8t8c8OjR\nG9+XJ3v+XOcPAE398pfSt751efm++/y35bgRnuz5AwCaunIOVDBzogh/ALhJTJlyeT6EZQX3UCBP\n39jNKy5ckI4fl7p391/uBQBuuO02/32Ptm2Tbr9dGjbsxvfFmH8LTp3y3+v8yBEpOVl64w1pwABX\nSwKAoDHs04JVq/zBL12e7QcANzvCvwUNDc0vA8DNiPBvQW6u1KuX/31CAj1/AO0DJ3xbkJoqHTjg\nn0zRq5cUF+d2RQAQPE+f8OUZvgDgDE+Hvxeu9gGA9ogxfwAwEOEPAAYi/AHAQIQ/ABiI8AcAAxH+\nAGAgwh8ADET4A4CBPB3+PMMXAJzBDF8AMJCne/4AAGcQ/gBgoKDDf8OGDRo7dqySkpLk8/m0b9++\nJtvcc8898vl8ja+wsDDNmTMn2KYBADco6PCvra3VyJEj9dxzz8myrKtuY1mWvvGNb6iiokLl5eU6\ndeqUVq5cGWzTAIAbFPTDXKZNmyZJKisrU3PnjiMjI5WUlBRscwCAEGizMf/f/e53SkpK0oABA/S9\n731P//3vf9uqaQDAFdrkMY5Tp05Vr1691L17d+3bt09PPPGEjhw5oj/96U9t0TwA4AqtCv/8/HzN\nnj1bkn8cf8uWLcrMzGzxc7NmzWp8379/f3Xr1k2jR4/Whx9+qLS0tFaWjLo66aWXpE8+kaZN8z9Y\nHgBao1Xhn5WVpREjRjQup6Sk3FCjw4cPl23bOnr0aLPhn52drfDwwBJ5nq+UlSVt3ep/n5cnlZRI\nUVHu1gTg5tKq8I+KilJ6evo1v36tq32utGfPHlmWpVtvvbXZ7QoKCpjhe4X//Ody8EvSoUPSnj3S\nyJHu1QTg5hP0mH9lZaWOHz+uf/3rX7JtW4cOHZJt2+rWrZuSk5N17Ngx5efna/z48erSpYv27t2r\nRYsW6e6779ZnP/vZUHwPRomL8w/zfPyxfzk8XOrRw92aANx8gr7aZ+PGjRo8eLAmTJggy7KUk5Oj\nIUOG6MUXX5Qk3XLLLXrzzTc1ZswYZWRk6Lvf/a4mT56sjRs3Bl28iSIipI0bpcGDpb59pd/+Vurd\n2+2qANxsuLEbABiIe/sAgIEIfwAwEOEPAAYi/AHAQJ4Ofx7jCADO4GofADCQp3v+AABnEP4AYCDC\nHwAMRPgDgIEIfwAwEOEPAAYi/AHAQIQ/ABjI0+HPDF8AcAYzfAHAQJ7u+QMAnEH4A4CBCH8AMBDh\nDwAGIvwBwECEPwAYiPAHAAMR/gBgIE+HPzN8AcAZzPAFAAN5uucPAHAG4Q8ABiL8AcBAhD8AGIjw\nBwADEf4AYCDCHwAMRPgDgIE8Hf7M8AUAZzDDFwAM5OmePwDAGUGFf319vRYvXqyBAwcqOjpaKSkp\nmjFjhk6dOhWwXWVlpaZOnaq4uDjFx8dr1qxZqq2tDapwAMCNCyr8z58/r9LSUj3zzDPas2ePNmzY\noMOHDysrKytguylTpujgwYMqLCzU5s2btX37ds2ePTuowgEANy7kY/4lJSUaPny4ysrK1KNHDx08\neFD9+/fX7t27NXjwYEnS1q1b9ZWvfEUnTpxQt27dmuyDMX8AcFbIx/yrqqpkWZY6d+4sSdq1a5fi\n4+Mbg1+SRo8eLcuyVFxcHOrmAQDXIaThX1dXpyVLlmjKlCmKjo6WJJWXl6tr164B24WFhSkhIUHl\n5eWhbB4AcJ1aFf75+fmKiYlRTEyMYmNjtWPHjsav1dfXa/LkybIsSy+88EKL+7JtW5Zltb5iAEDQ\nwluzcVZWlkaMGNG4nJKSIuly8H/00Ud66623Gnv9ktStWzedPn06YD8XL15UZWWlkpOTm20vOztb\n4eGBJebk5CgnJ6c1ZQMArhD0Cd9LwX/s2DEVFRUpISEh4OuHDh1S//79VVJS0jju//rrr2v8+PGc\n8AUAlwQV/hcvXtTEiRNVWlqq1157LWBsPyEhQREREZKk8ePH6/Tp08rLy9OFCxf09a9/XcOGDdPL\nL7981f0S/gDgrKDCv6ysTOnp6QHrLo3lFxUV6a677pLkvwJo3rx52rRpk3w+nyZNmqSf/exnioyM\nvOp+CX8AcBb39gEAA3FvHwAwEOEPAAYi/AHAQIQ/ABiI8AcAA3k6/HmMIwA4g0s9AcBAnu75AwCc\nQfgDgIEIfwAwEOEPAAYi/AHAQIQ/ABiI8AcAAxH+AGAgT4c/M3wBwBnM8AUAA3m65w8AcAbhDwAG\nIvwBwECEPwAYiPAHAAMR/gBgIMIfAAxE+AOAgTwd/szwBQBnMMMXAAzk6Z4/AMAZhD8AGIjwBwAD\nEf4AYCDCHwAMRPgDgIEIfwAwkKfDn0leAOAMJnkBgIE83fMHADiD8AcAAwUV/vX19Vq8eLEGDhyo\n6OhopaSkaMaMGTp16lTAdr1795bP52t8hYWFaeXKlUEVDgC4ceHBfPj8+fMqLS3VM888o4EDB6qy\nslLz589XVlaW3n333cbtLMvS0qVL9fjjj+vSKYaYmJjgKgcA3LCgwj82NlZbt24NWPeLX/xCw4cP\n14kTJ9SjR4/G9dHR0UpKSgqmOQBAiIR8zL+qqkqWZalz584B61esWKHExEQNGTJEq1at0sWLF0Pd\nNADgOgXV879SXV2dlixZoilTpig6Orpx/YIFCzRkyBAlJCTonXfe0ZIlS1ReXq5Vq1aFsnkAwHVq\n1XX++fn5mj17tv+DlqUtW7YoMzNTkv/k78SJE3Xq1CkVFRUFhP+VfvOb3+ib3/ymzp07p4iIiCZf\n5zp/AHBWq8K/trZWFRUVjcspKSnq0KGD6uvrNXnyZP3zn//UW2+9pfj4+Gb3c+DAAQ0YMECHDh3S\n7bff3uTrl8J/3LhxCg8P/OMkJydHOTk511syAOAqWjXsExUVpfT09IB1l4L/2LFjKioqajH4JWnP\nnj3y+Xzq2rVrs9sVFBTQ8wcABwQ15n/x4kU9/PDDKi0t1WuvvaZPP/208S+DhIQERUREaNeuXSou\nLtaoUaMUExOjd955R4sWLdL06dMVFxcXkm8CANA6Qd3bp6ysrMlfArZty7IsFRUV6a677tKePXs0\nZ84cHT58WHV1dUpLS9PXvvY15ebmXnW8X2LMHwCcxo3dAMBA3NsHAAxE+AOAgQh/ADAQ4Q8ABiL8\nAcBAng5/nuELAM7gUk8AMJCne/4AAGcQ/gBgIMIfAAxE+AOAgTx5wte2bZ09e1YxMTGyLMvtcgCg\n3fFk+AMAnMWwDwAYiPAHAAMR/gBgIMIfAAxE+AOAgQh/ADAQ4Q8ABvo/8m+fkcHrK6EAAAAASUVO\nRK5CYII=\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# 1/β_ML = 1/N Σ{y(x_n,W_ML) - t_n}^2から計算したβ_MLでエビデンスを計算\n", "beta_ML_List = []\n", "ln_p_List = []\n", "for M in range(10):\n", " Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1))] for x in X.list()]) \n", " Phi_t = Phi.transpose()\n", " # m_N, Aを定義\n", " A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi\n", " m_N = beta* A.inverse() * Phi_t* t \n", " # エビデンスの対数\n", " res = (t - Phi*m_N)\n", " beta_ML = N / (res*res)\n", " E_mN = beta_ML/2* res * res + alpha/2 * m_N * m_N\n", " ln_p_t = M/2 * ln(alpha) + N/2 * ln(beta_ML) - E_mN - 1/2 * ln(A.det()) - N/2 * ln(2*pi)\n", " ln_p_List += [ln_p_t]\n", " beta_ML_List += [beta_ML]\n", "list_plot(ln_p_List, ymin = -26, ymax = 0).show(figsize=4)\n", "#print beta_ML_List" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t\n", "\t\tどうして図3.14のような結果にならないのかと調べたのですが、他にも同様の計算をした\n", "\t\t方がいらして、私と同じ傾向になったとの記述がありました。\n", "\t\t基底関数を色々変えて、線形回帰のエビデンスを計算してみた\n", "\t
\n", "\t\n", "\t\t私の推測では、ln |A|の計算をA.norm()と間違えたのではないかと思われます。以下にA.norm()としたときの図を示します。\n", "\t
\n", "" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEdCAYAAADkeGc2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAHUNJREFUeJzt3X1cVvX9x/H3wYtMuRPSQLFUsiVa4E1TEmdWNm9KWTqM\nSzOr2SinlralW030Uc1pd2s/tfpZa66FuM1YhPPXL4XKH3NO1Au7EaoxMZU0HwM1cMTN+f3BxAhE\n4bouzpHzej4ee8Q517m+3w/Xo7378jlfDoZpmqYAAI4SYHUBAID2R/gDgAMR/gDgQIQ/ADgQ4Q8A\nDkT4A4ADEf4A4ECEPwA4EOEPAA5E+AOAA9ky/FevXq1+/fqpS5cuSkhI0M6dO60uyTaWL1+u4cOH\nKzQ0VJGRkbrtttv08ccfW12WbS1fvlwBAQFauHCh1aXYzuHDhzVz5kx1795dXbt2VXx8vHbv3m11\nWbZRV1enn//854qJiVHXrl3Vv39/Pf7441aX5TO2C/8NGzbooYce0rJly7Rnzx7Fx8dr3LhxOnbs\nmNWl2cK2bds0b9487dixQ1u2bFF1dbW++93v6tSpU1aXZjs7d+7U2rVrFR8fb3UptlNeXq7ExER1\n7txZb731lvbt26enn35a4eHhVpdmG7/85S/14osvas2aNSosLNTKlSu1cuVKrVq1yurSfMO0mREj\nRpjz589vOK6rqzOjo6PNFStWWFiVfX3xxRemYRjmtm3brC7FVk6ePGl+61vfMrdu3WqOGTPGXLBg\ngdUl2cqiRYvM0aNHW12Grd16663m7NmzG52bOnWqOXPmTIsq8i1brfyrq6u1a9cu3XTTTQ3nDMPQ\n2LFjtX37dgsrs6/y8nIZhqGIiAirS7GVH/3oR5o0aZJuvPFGq0uxpTfffFPXXnutpk2bpsjISA0d\nOlQvvfSS1WXZysiRI7V161Z98sknkqSCggLl5eVp4sSJFlfmGy6rC/i6Y8eOqba2VpGRkY3OR0ZG\nqqioyKKq7Ms0TT344IMaNWqUBg4caHU5tpGRkSGPx6P8/HyrS7Gt4uJiPf/883rooYf0yCOPaMeO\nHZo/f74uvvhi3XHHHVaXZwuLFy/WiRMnNGDAAHXq1El1dXV64oknlJKSYnVpPmGr8D8b0zRlGIbV\nZdjOnDlz9NFHHykvL8/qUmzj4MGDevDBB/X2228rMDDQ6nJsq66uTsOHD9djjz0mSYqPj9eHH36o\n559/nvD/jw0bNig9PV0ZGRkaOHCgPB6PHnjgAfXq1UszZ860ujyv2Sr8u3fvrk6dOunIkSONzh89\nerTJTwNON3fuXP3lL3/Rtm3b1LNnT6vLsY1du3bpiy++0LBhw2T+5+8U1dbW6r333tOqVatUVVXF\nQkJSz549FRsb2+hcbGysXn/9dYsqsp+HH35YP/vZz5ScnCxJGjRokPbv36/ly5d3iPC3Vc8/MDBQ\nw4YN09atWxvOmaaprVu3auTIkRZWZi9z587VG2+8odzcXF1++eVWl2MrY8eO1fvvvy+Px6OCggIV\nFBTo2muv1R133KGCggKC/z8SExObtFKLiorUp08fiyqyn8rKyib/vgQEBKiurs6iinzLVit/SVq4\ncKFmzZqlYcOGafjw4Xr22WdVWVmpu+66y+rSbGHOnDlav369srKyFBQU1PBTUlhYmC6++GKLq7Ne\nUFBQk/sfQUFBuuSSS5qsdJ1swYIFSkxM1PLlyzVt2jTt2LFDL730ktauXWt1abYxadIkPfHEE7rs\nsss0aNAg7d69W88++6xmz55tdWm+Ye1mo+atXr3a7NOnj3nxxRebCQkJ5s6dO60uyTYMwzADAgKa\n/G/dunVWl2ZbN9xwA1s9m7Fp0ybzmmuuMbt06WIOHDjQfPnll60uyVa+/PJLc8GCBWbfvn3Nrl27\nmv379zeXLFliVldXW12aTximyR9wBwCnsVXPHwDQPgh/AHAgwh8AHIjwBwAHIvwBwIEIfwBwIEt/\nycs0TZ08edLKEgCgwwkJCTnnb7NbGv4nT55UWFiYlSUAQIdz/PhxhYaGtniNpb/k1dLKPyUlRRkZ\nGV7P4atxfDlWRx3Hl2N11HF8OVZHHceXY3XUcc41lu1X/oZhnPW/Ti6X65z/5TofvhrHl2N11HF8\nOVZHHceXY3XUcXw5Vkcdxxdjdfgbvm6323Zj+bImX+AzOjc+o3PjMzo3O31Gtn22z+TJk5WVlWV1\nGbbGZ3RufEbnxmd0bh3xM+rwK38AQFOdli5dutTqIs7mmmuusboE2+MzOjc+o3PjMzq3jvYZ2bbt\nAwDwH9o+AOBAhD8AOBDhDwAORPgDgAMR/gDgQLYMf9M0deLECbERCQD8w5bhf/ppnzzuGQD8w5bh\nDwDwL8IfPrFqlTRmjHTvvdKJE1ZXA+BcLH2k87mkpKTI5XLJ7Xbb7ul8OOPPf5bmzav/+t13pcpK\n6bXXrK0JQMtsHf4ZGRk+e/Y1/KegoPHx3r3W1AHg/Hnd9snMzNT48ePVo0cPBQQEaO83/p9fVlam\n+fPna8CAAQoKClKfPn30wAMP6AS9gQ7j5pulTp3OHI8fb10tAM6P1yv/iooKjRo1StOmTdO9997b\n5PXDhw+rtLRUzzzzjGJjY1VSUqLU1FSVlpbqD3/4g7fTwwZGjpTeflt64w3pyiul+++3uiIA5+Kz\np3qWlJSoX79+8ng8iouLa/HaP/3pT5o5c6YqKioUEND0h48TJ04oLCzsvP4IMQCg9SzZ7VNeXq7Q\n0NBmgx8A4H/tnr7Hjh3T448/rtTU1PaeGgDwH60K//T0dIWEhCgkJEShoaHKy8tr1WQnT57ULbfc\noquvvlppaWmtei8AwHdadcM3KSlJCQkJDcfR0dHn/d4vv/xS48aNU7du3fT666+r09e3h5zF6X3+\nX8eefwDwXqvCPygoSDExMWd93TCMZs+fPHlS48aNU5cuXZSVlaWLLrrovOZjnz8A+IfXWz3Lysp0\n4MABHTp0SKZpqrCwUKZpKioqSpGRkfryyy91880369///rdee+01lZeXN7z39O8GAADal9dbPdet\nW6e77767yao/LS1NS5Ys0bvvvqsbb7yx0WumacowDP3zn//U5Zdf3mRMtnoCgH/5bJ+/LxH+AOBf\n9FwAwIEIfwBwIFuHf0pKiiZPnqz169dbXQoAdCj0/AHAgWy98gcA+AfhDwAORPgDgAMR/gDgQIQ/\nADiQrf+AO9BWp05J//VfUlmZdNdd0lVXWV0RYC+23uo5YcIEuVwuHuOMVrv1VmnTpvqvIyKk99+X\nevWytibATmwd/uzzR1vU1UmBgfX/PG3jRmnKFOtqAuyGnj86nIAAKTb2zLHLRdsH+CbCHx1SVpY0\nebL0ne9IGRnSoEFWVwTYC20fAHAgVv4A4ECEPwA4kK3Dn0c6A4B/0PMHAAey9cofAOAfXod/Zmam\nxo8frx49eiggIEB79+5t8foJEyYoICBAWVlZ3k4NAGgjr8O/oqJCo0aN0ooVK2QYRovXPvvss+rU\nqdM5rwMA+JfXD3a74447JEklJSVq6fZBQUGBfvWrX2nnzp2KiorydloAgBfaped/6tQpTZ8+XatX\nr9all17aHlMCAFrQLuG/YMECjRo1Srfeemt7TAcAOIdWtX3S09OVmpoqSTIMQ5s3b1ZiYmKL78nK\nylJOTo48Hk+ri0tJSZHL1bjE9n68c12ddP/90p/+JF1xRf1zYmJi2m16APCLVoV/UlKSEhISGo6j\no6PP+Z7c3FwVFxcrLCys0fkpU6Zo9OjRysnJOet7MzIyLN/nv26d9N//Xf/1v/4lpaZKb79taUkA\n4LVWhX9QUJBiWlj2NreL56c//anuvffeRueuvvpqPffccxdEG6i0tOVjALgQeb3bp6ysTAcOHNCh\nQ4dkmqYKCwtlmqaioqIUGRmpSy+9tNmbvJdddpn69Onj7fR+l5wsPfVU/Z8DlOpX/gBwofP6hm9W\nVpaGDBmiSZMmyTAMud1uDR06VC+++OJZ33Mh7fO/8kppzx7p5Zeld96R5s2zuiIA8B7P9gEAB+LZ\nPgDgQIQ/ADiQrcOf5/kDgH/Q8wcAB7L1yh8A4B+EPwA4EOEPAA5E+AOAAxH+AOBAtg5/tnoCgH+w\n1RMAHMjWK38AgH8Q/gDgQIQ/ADgQ4Q8ADkT4A4ADEf4A4EC2Dn/2+aMjOXVK2r1bOnrU6koA9vkD\n7eLYMek735EKC6WuXaU//1m6+Warq4KTeb3yz8zM1Pjx49WjRw8FBARo7969zV63fft23XTTTQoO\nDlZYWJjGjBmjqqoqb6cHLghr19YHvyRVVkqPPmptPYDX4V9RUaFRo0ZpxYoVMgyj2Wu2b9+uCRMm\naPz48crPz1d+fr7mzp2rgABbd50An3G5Gh8HBlpTB3Caz9o+JSUl6tevnzwej+Li4hq9dt1112nc\nuHFaunTpeY1F2wcdzYkT0tix0s6dUkSElJ0tXXed1VXByfy+9P7iiy+0Y8cOde/eXYmJiYqKitKY\nMWOUl5fn76kB2wgNlbZvl0pKpM8+I/hhPb+Hf3FxsSRp2bJlSk1N1VtvvaWhQ4fqpptu0j/+8Q9/\nTw/YRqdO0uWX19/wBazWqvBPT09XSEiIQkJCFBoael6r97q6OknSfffdpzvvvFPx8fF65plndNVV\nV+k3v/lN26oGAHjFde5LzkhKSlJCQkLDcXR09Dnf07NnT0lSbGxso/OxsbE6cOBAi+9NSUmR6xt3\nytxut9xu9/mWDABoRqvCPygoSDExMWd9vbndPn379lWvXr1UVFTU6PzHH3+siRMntjhfRkYGN3wB\nwA9aFf7NKSsr04EDB3To0CGZpqnCwkKZpqmoqChFRkZKkn7yk59o6dKliouL0+DBg/Xb3/5WRUVF\n2rhxo9ffAACg9bwO/6ysLN19990yDEOGYTS0ZNLS0rRkyRJJ0gMPPKCqqiotXLhQ//rXvxQfH68t\nW7aoX79+3k4PAGgDHu8AAA7Er9gCgAMR/gDgQLYOfx7pfHb/939SVpZUUWF1JQAuRPT8L0CPPCL9\n4hf1X8fHS3l5UlCQtTUBuLDYeuWPpkxTevrpM8cFBdLbb1tXD4ALE+F/gTEMKSys8bnwcGtqAXDh\nIvwvQK+9JvXoUf9M+J/8RLr+eqsrAnChoed/Aaurk/h7OADagui4gBH8ANqK+AAAB7J1+LPPHwD8\ng54/ADiQrVf+AAD/IPwBwIEIfwBwIMIfAByI8AcAB7J1+LPVEwD8g62eAOBAXq/8MzMzNX78ePXo\n0UMBAQHau3dvk2uOHDmimTNnqmfPngoODtawYcP0+uuvezs1AKCNvA7/iooKjRo1SitWrJBhGM1e\nM3PmTH3yySfKzs7WBx98oClTpmjatGkqKCjwdnoAQBv4rO1TUlKifv36yePxKC4urtFrISEheuGF\nFzRjxoyGc927d9fKlSt1zz33NBmLtg8A+Fe73PBNTEzUhg0bVFZWJtM0lZGRoaqqKo0ZM6Y9pgcA\nfIOrPSbZsGGDbr/9dl1yySVyuVwKCgpSZmamYmJi2mN6AMA3tGrln56erpCQEIWEhCg0NFR5eXnn\n9b5HH31Ux48fV05Ojnbt2qWFCxcqOTlZH374YZuKBgB4p1U9/4qKCh05cqThODo6Wp07d5Z09p5/\ncXGx+vfvr48++kgDBgxoOH/zzTfryiuv1Jo1a5rMc7rnP2HCBLlcjX84cbvdcrvd5/8dAgCaaFXb\nJygoqMVWTXO7fSorK2UYRpPXOnXqpLq6uhbny8jI4IYvAPiB1z3/srIyHThwQIcOHZJpmiosLJRp\nmoqKilJkZKQGDBigK664QqmpqXryySd1ySWXKDMzU1u2bNGmTZt88T0AAFrJ690+WVlZGjJkiCZN\nmiTDMOR2uzV06FC9+OKLkiSXy6XNmzerR48emjx5suLj4/X73/9ev/vd7zRu3DivvwEAbXP8uPTI\nI9KcOVIzv5uJDo7HOwAOdcMN0jvv1H8dFiZ98IHUu7elJaEd2frBbgD8o6bmTPBL9T8F7NplWTmw\nAOEPOJDLJcXHnznu3FkaNMi6etD+CH/AobKzpRkzpFtukd58U+rf3+qK0J5s3fM/vc+fvf0A4Fu2\nDn9u+AKAf9D2AQAHIvwBwIEIfwBwIMIfAByI8AcAB7J1+KekpGjy5Mlav3691aUAQIfCVk8AcCBb\nr/wBAP5B+AOAAxH+AOBAhD8AOBDhDwAORPgDgAPZOvzZ5w8A/sE+fwBwIK9W/jU1NVq0aJHi4uIU\nHBys6OhozZo1S6WlpY2uKysr04wZMxQWFqbw8HDNnj1bFRUVXhUOAGg7r8K/srJSHo9HaWlp2rNn\njzIzM1VUVKSkpKRG102fPl379u3T1q1btWnTJr333ntKTU31qnAAQNv5vO2Tn5+vESNGqKSkRL17\n99a+ffs0aNAg7dq1S0OGDJEkvfXWW7rlllt08OBBRUVFNRmDtg8A+JfPb/iWl5fLMAx169ZNkvS3\nv/1N4eHhDcEvSWPHjpVhGNqxY4evpwcAnAefhn9VVZUWL16s6dOnKzg4WJL0+eef69JLL210XadO\nnRQREaHPP//cl9MDAM6TqzUXp6enN/TqDcPQ5s2blZiYKKn+5m9ycrIMw9CaNWvOOZZpmjIMo8Vr\nUlJS5HI1LtHtdsvtdrembADAN7Qq/JOSkpSQkNBwHB0dLelM8H/22WfKyclpWPVLUlRUlI4ePdpo\nnNraWpWVlSkyMrLF+TIyMuj5A4AftCr8g4KCFBMT0+jc6eAvLi5Wbm6uwsPDG71+3XXXqby8XHv2\n7Gno+2/dulWmaWrEiBFelg8AaAuvdvvU1tZqypQp8ng8ys7ObtTbj4iIUGBgoCRp4sSJOnr0qJ5/\n/nl99dVXuueeezR8+HC9+uqrzY7Lbh8A8C+vwr+kpKTJTwKne/m5ubkaPXq0pPodQHPnztWbb76p\ngIAAff/739dzzz2nrl27Njsu4Q8A/mW7xzsUFko//OEJbdsWpl//+rjmzSP8AcDXbBf+AwdK+/ad\nkBQm6bjy80M1bJjVVQFAx2K7p3r+4x8tHwMAvGe78P/+9898fdFFKXr5ZR7pDAC+Zru2T02NtHr1\nCT34YJj27j2ua66h5w8Avma78JfY7QMA/ma7tg8AwP8IfwBwIMIfgG1kZUk/+IH05JP19//gP616\ntg8A+MuWLdL3viedvgv5+efS009bW1NHxsofgC28886Z4Jek3FzLSnEEW4d/SkqKJk9mnz/gBNde\n2/IxfIutngBsY+1a6Y03pKuukh5/XOrSxeqKOi7CHwAcyNZtHwCAfxD+AOBAhD8AOBDhDwAOZOvw\nZ6snAPgHu30AwIG8WvnX1NRo0aJFiouLU3BwsKKjozVr1iyVlpY2XFNSUqLZs2crJiZGXbt21ZVX\nXqmlS5equrra6+IBAG3j1bN9Kisr5fF4lJaWpri4OJWVlWn+/PlKSkrS3//+d0lSYWGhTNPU2rVr\ndcUVV+iDDz7Q7NmzVVlZqZUrV/rkmwAAtI7P2z75+fkaMWKESkpK1Lt372aveeqpp/TCCy/o008/\nbfZ12j4A4F8+v+FbXl4uwzDUrVu3Fq+JiIjw9dQAgPPk0/CvqqrS4sWLNX36dAUHBzd7zaeffqpV\nq1bpvvvu8+XUAIBWaFX4p6enKyQkRCEhIQoNDVVeXl7DazU1NUpOTpZhGFqzZk2z7z906JAmTJig\n22+/Xffcc493lQMA2qxVPf+KigodOXKk4Tg6OlqdO3duCP79+/crJydH4eHhTd57+PBh3XDDDRo5\ncqReeeWVFuc53fOfMGGCXK7G96Tdbrfcbvf5lgwAaIbXN3xPB39xcbFyc3Ob7eUfOnRIN954o779\n7W/r1VdflWEYLY7JDV8A8C+vtnrW1tZq6tSp8ng8ys7OVnV1dcNPBhEREQoMDFRpaanGjBmjvn37\nauXKlTp69GjD+yMjI72rHgDQJl6t/EtKShQTE9PonGmaMgxDubm5Gj16tNatW9ekv3/6mtra2mbH\nZeUPAP7F4x0AwIFs/WA3AIB/EP4A4EC2Dn8e6QwA/kHPHwAcyNYrfwCAfxD+AOBAhD8AOBDhDwAO\nRPgDgAMR/gDgQLYOf/b5A4B/sM8fABzI1it/AIB/EP4A4ECEPwA4EOEPAA5E+AOAAxH+AHAWJSVS\nSoo0frz0v/9rdTW+ZeutnhMmTJDL5ZLb7Zbb7ba6LAAOc/XV0ocf1n/dubP0wQdS//7W1uQrtg5/\n9vkDsEp1tXTRRY3PZWVJkyZZU4+vedX2qamp0aJFixQXF6fg4GBFR0dr1qxZKi0tbfb6r776SoMH\nD1ZAQID27t3rzdQA4FeBgdL11585Dg+Xhg2zrh5f8yr8Kysr5fF4lJaWpj179igzM1NFRUVKSkpq\n9vqHH35YvXv3lmEY3kwLAO0iK0t69FFp3jxp2zapVy+rK/Idn7d98vPzNWLECJWUlKh3794N5zdv\n3qwf//jH2rhxowYOHCiPx6O4uLhmx6DtAwD+5fL1gOXl5TIMQ926dWs4d+TIEf3whz9UVlaWunTp\n4uspAQCt5NOtnlVVVVq8eLGmT5+u4ODghvN333235syZoyFDhvhyOgBAG7Vq5Z+enq7U1FRJkmEY\n2rx5sxITEyXV3/xNTk6WYRhas2ZNw3t+/etf6+TJk1q0aJEkqTVdppSUFLlcjUtk2ycAeK9VPf+K\nigodOXKk4Tg6OlqdO3duCP79+/crJydH4eHhDdfcdtttys7ObjRObW2tXC6XZsyYoVdeeaXJPPT8\nAcC/vL7hezr4i4uLlZubq4iIiEavHzx4UCdOnGg4Pnz4sMaNG6eNGzdq+PDh6tXM7XPCHwD8y6sb\nvrW1tZo6dao8Ho+ys7NVXV3d8JNBRESEAgMDG+34kaSgoCCZpqmYmJhmgx8A4H9ehf/BgwcbWjqD\nBw+WVN/TNwxDubm5Gj16dLPvY58/AFiLxzsAgAPxVE8AcCDCHwAcyNbhn5KSosmTJ2v9+vVWlwIA\nHQo9fwBwIFuv/AEA/kH4A4ADEf4A4ECEPwA4EOEPAA5E+AOAA9k6/NnnDwD+wT5/AHAgW6/8AQD+\nQfgDgAMR/gDgQIQ/ADgQ4Q8ADmTr8GerJwD4B1s9AcCBvFr519TUaNGiRYqLi1NwcLCio6M1a9Ys\nlZaWNrl206ZNSkhIUNeuXRUREaEpU6Z4MzUAwAtehX9lZaU8Ho/S0tK0Z88eZWZmqqioSElJSY2u\n27hxo+6880794Ac/0Pvvv6+//vWvmj59uleFAwDazudtn/z8fI0YMUIlJSXq3bu3amtr1bdvXz32\n2GO66667zmsM2j4A4F8+v+FbXl4uwzDUrVs3SdLu3bt1+PBhSdLQoUPVq1cvTZw4UR999JGvpwYA\nnCefhn9VVZUWL16s6dOnKzg4WJJUXFws0zS1bNkyLVmyRJs2bVJ4eLiuv/56lZeX+3J6AMB5alX4\np6enKyQkRCEhIQoNDVVeXl7DazU1NUpOTpZhGFqzZk3D+bq6OknSo48+qu9973saMmSIXnnlFRmG\noT/+8Y8++jYAAK3has3FSUlJSkhIaDiOjo6WdCb4P/vsM+Xk5DSs+iWpZ8+ekqTY2NiGcxdddJFi\nYmJ04MCBFudLSUmRy9W4RLfbLbfb3ZqyAQDf0KrwDwoKUkxMTKNzp4O/uLhYubm5Cg8Pb/T6sGHD\n1LlzZxUVFWnkyJGSpOrqau3fv199+vRpcb6MjAxu+ALA12RmSllZUmystHCh5GpVip/RxrfVq62t\n1dSpU+XxeJSdna3q6modOXJEkhQREaHAwECFhITovvvuU1pamnr37q0+ffpo5cqVMgxDycnJ3kwP\nAI7yP/8jff1XpI4elZ56qm1jeRX+Bw8eVHZ2tiRp8ODBkiTTNGUYhnJzczV69GhJ0lNPPaXAwEDd\neeedOnXqlEaMGKGcnByFhYV5Mz0AOMp77zU+fvfdto/F4x0A4AKRmdl45X///dLX9te0ilcrfwBA\n+7ntNumll6Q33qjv+S9b1vaxWPkDgAPxSGcAcCBW/gDgQLYMf9M0dfLkSYWEhMgwDKvLAYAOx5bh\nDwDwL1v3/AEA/kH4A4ADEf4A4ECEPwA4EOEPAA5E+AOAAxH+AOBA/w+ACIOuvTjNCAAAAABJRU5E\nrkJggg==\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 当初、式3.85, 3.86のミスプリと考えたのですが、\n", "# →http://d.hatena.ne.jp/n_shuyo/20090715/evidence#c\n", "# n_shuyoさんのコメントで 正則行列の行列式は、逆行列の行列式の逆数:|A| = 1/|A^-1|\n", "# 式3.85, 3.86 は正しい!\n", "# 私の推測:図3.14の計算で ln |A|とするところをln A.norm()と間違えたのではないか\n", "ln_p_List = []\n", "for M in range(10):\n", " Phi = matrix(RDF, [[ _phi(x,j) for j in range(0, (M+1))] for x in X.list()]) \n", " Phi_t = Phi.transpose()\n", " # m_N, Aを定義\n", " A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi\n", " m_N = beta* A.inverse() * Phi_t* t \n", " # エビデンスの対数\n", " res = (t - Phi*m_N)\n", " E_mN = beta/2* res * res + alpha/2 * m_N * m_N\n", " ln_p_t = M/2 * ln(alpha) + N/2 * ln(beta) - E_mN - 1/2 * ln(A.norm()) - N/2 * ln(2*pi)\n", " ln_p_List += [ln_p_t]\n", "list_plot(ln_p_List, figsize=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t\n", "\t\tα、βを固定にした場合に、エビデンスが最大はM=4だったので、M=4に対する\n", "\t\t最適なα、βを求めてみます。\n", "\t
\n", "\t\n", "\t\tエビデンスを最大にするαは、式(3.92)\n", "$$\n", "\t\t\\alpha = \\frac{\\gamma}{m_N^T m_N}\n", "$$\n", "\t\tここで$\\gamma$は、式(3.87)\n", "$$\n", "\t\t(\\beta \\Phi^T \\Phi) u_i = \\lambda_i u_i\n", "$$\n", "\t\tの固有値から、式(3.91)\n", "$$\n", "\t\t\\gamma = \\sum_{i} \\frac{\\lambda_i}{\\alpha + \\lambda_i}\n", "$$\n", "\t\tで求まります。\n", "\t
\n", "\t\n", "\t\tβについても、式(3.95)\n", "$$\n", "\t\t\\frac{1}{\\beta} = \\frac{1}{N - \\gamma} \\sum_{n=1}^{N} ( t_n - m_N^T \\phi(x_n))^2\n", "$$\n", "\t
\n", "\t\n", "\t\t$\\gamma$は、$\\alpha$に依存し、$m_N$も$\\alpha$に依存するため、すぐに最適な値を得ることが\n", "\t\tできません。\n", "\t\t\n", "\t\t
\n", "\t\t以下の例では、\n", "\t\t