![](\"images/03/Figure3.8a.png\")
![](\"images/03/Figure3.8b.png\")
![](\"images/03/Figure3.8c.png\")
![](\"images/03/Figure3.8d.png\")
\n", "\t\t第3章-Sageを使って逐次ベイズ学習を試してみるでは、\n", "\t\t重みwの分布について解いていましたが、\n", "\t\t実際には新しいｘの値たいするtを予測することが目的となります。\n", "\t\tここではPRML3.3.2の予測分布、図3.8(PRMLのサイトから引用)をSageを使って試してみます。\n", "\t\t\n", "\t\t
| \n",
"\t\t\t\t | \n",
"\t\t\t
| \n",
"\t\t\t\t | \n",
"\t\t\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.8の再現\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", "x = var('x')\n", "sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green')\n", "data_plt = list_plot(zip(X, t))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t\n", "\t\t第3章-エビデンス近似をSageで試す\n", "\t\tのようにガウス基底関数を定義しますが、今回は0から1の範囲で定義します。\n", "\t
\n", "\t\n", "\t\t近似に使う場合には、_phiのようにj=0が1となるような項を追加します。\n", "\t
\n", "" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from pylab import linspace\n", "# 定数をセット\n", "M=9\n", "mu = linspace(0, 1, M)\n", "s = 0.1\n", "s_sq = (s)^2\n", "# ガウス基底関数\n", "def _phi_gauss(x, j):\n", " return e^(-1*(x - mu[j])^2/(2* s_sq))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# ガウス基底関数で三角関数の例題を近似\n", "# j=0に対応するため_phiを以下のように定義\n", "def _phi(x, j):\n", " if j == 0:\n", " return 1\n", " else:\n", " return _phi_gauss(x, j-1)\n", "# 初期化\n", "alpha = 2\n", "beta = 25" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t\n", "\t\t予測分布は、式(3.58)\n", "$$\n", "\t\tp(t|x, t, \\alpha, \\beta) = \\mathcal{N}(t|m_N^T\\phi(x), \\sigma_N^2(x))\n", "$$\n", "\t\tから計算し、その分散$\\sigma_N^2(x)$は、\n", "$$\n", "\t\t\\sigma_N^2(x) = \\frac{1}{\\beta} + \\phi(x)^T S_N \\phi(x)\n", "$$\n", "\t\tで与えられます。\n", "\t
\n", "" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 予測分布をプロットする\n", "def _plot_predict(X, t):\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", " Phi_dag = (alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi).inverse()*Phi_t;\n", " # 平均の重み\n", " Wml = beta*Phi_dag * t\n", " f = lambda x : (sum((Wml[i]*_phi(x, i)) for i in range(0, (M+1))))\n", " # 分散(標準偏差)\n", " def s(x):\n", " phi_x = vector([_phi(x, i).n() for i in range(M+1)])\n", " S = (alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi).inverse()\n", " s_sqr = 1/beta + phi_x * S * phi_x\n", " return sqrt(s_sqr)\n", " data_plt = list_plot(zip(X, t))\n", " s_u_plt = plot(lambda x : f(x) + s(x), [x, 0, 1], rgbcolor='grey')\n", " s_d_plt = plot(lambda x : f(x) - s(x), [x, 0, 1], rgbcolor='grey')\n", " y_plt = plot(f, [x, 0, 1], rgbcolor='red')\n", " (y_plt + data_plt + sin_plt + s_u_plt + s_d_plt).show(ymin=-1.5, ymax=1.5, figsize=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t\n", "\t\t以下にM=1, 2, 4, 10の予測分布を示します。パラメータは、\n", "\t\t