{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "オリジナル作成:　2011/06/04" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\tHiroshi TAKEMOTO\n", "\t(take.pwave@gmail.com)\n", "\t\n", "\t

# 第3章-エビデンス近似をSageで試す

\n", "\t

\n", "\t\t第1章-Sageを使って線形回帰を試してみるでは、\n", "\t\t与えられたモデルパラメータM, α、βに対して解いていましたが、\n", "\t\t実際には自分でベストなモデルパラメータM, α、βを求める必要があります。\n", "\t\t\n", "\t\tパターン認識と機械学習\n", "\t\tの３章ではエビデンス関数を使って最適なパラメータを求める方法が説明されています。\n", "\t

\t\n", "\t

\n", "\t\tここでは、Sageを使ってエビデンス関数の評価、最適なモデルパラメータの推定を試してみます。\n", "\t

\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

## エビデンス関数の評価

\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

\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

\n", "\t\tまずは、α、βを固定値\n", "\t\t

\n", "\t\t\t
• $\\alpha = 5 \\times 10^{-3}$
• \n", "\t\t\t
• $\\beta = 11.1$
• \n", "\t\t
\n", "\t\tとして、多項式の次数をMを0から9まで変えた値を計算し、プロットしてみます。\n", "\t

\n", "\t

\n", "\t\t一番良いエビデンスは、M=4あたりになります。\n", "\t\t理由は分かりませんが、図3.14のようなM=3以降の急激な下降は見られません。\n", "\t

### β_MLを使ってエビデンスを計算

\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

### どうして図3.14と合わないのか？

\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", "\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\t
1. 適当な $\\alpha, \\beta$を初期値とする
2. \n", "\t\t\t
3. $m_N$, $\\gamma$を求め
4. \n", "\t\t\t
5. 式(3.95)から$\\beta$を求める
6. \n", "\t\t\t
7. 式(3.92)から$\\alpha$を求め、2番目からを繰り返す
8. \n", "\t\t
\n", "\t

\n", "\t

\n", "\t\t以下の例では、\n", "\t\t

\n", "\t\t\t
• $\\alpha = 5 \\times 10^{-3}$
• \n", "\t\t\t
• $\\beta = 11.1$
• \n", "\t\t
\n", "\t\tから20回の計算を繰り返し、$\\gamma, \\alpha, \\beta$の値をプリントアウトしてみました。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4.02295319874 0.0160029572678 15.5703330431\n", "3.86957616939 0.0192233025727 14.6079344569\n", "3.81804019954 0.0208658867531 14.0162355337\n", "3.79079625631 0.021876053777 13.6574090928\n", "3.77410819657 0.0225480023858 13.4248103487\n", "3.76308615077 0.0230152846069 13.2666278994\n", "3.75546956914 0.0233495381967 13.1554382026\n", "3.75004882783 0.0235932327887 13.0754487343\n", "3.74611222743 0.0237732924978 13.0169432556\n", "3.74321234677 0.0239076177514 12.973632867\n", "3.7410540217 0.0240085312871 12.9412853436\n", "3.73943541844 0.0240847392383 12.9169657574\n", "3.7382147331 0.0241425141272 12.8985911035\n", "3.73729026269 0.0241864427411 12.8846563016\n", "3.73658790674 0.0242199173426 12.8740587014\n", "3.73605302024 0.0242454685376 12.8659818046\n", "3.73564493065 0.0242649966702 12.8598160009\n", "3.73533314922 0.0242799360403 12.8551032449\n", "3.73509469602 0.0242913734086 12.8514976848\n", "3.73491217789 0.0243001346625 12.848737196\n" ] } ], "source": [ "# M=4の場合のα、βを決める\n", "M=4\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", "# Φ^T Φは固定なので、先に固有値を計算\n", "B = Phi_t*Phi\n", "lambs = B.eigenvalues()\n", "# 初期化\n", "alpha = 5*10^(-3)\n", "beta = 11.1\n", "gamma = M\n", "for i in range(20):\n", " A = alpha*matrix((M+1),(M+1),1) + beta*Phi_t * Phi\n", " m_N = beta* A.inverse() * Phi_t* t \n", " gamma = sum((lamb*beta) /(alpha + (lamb*beta)) for lamb in lambs)\n", " alpha = gamma/(m_N*m_N)\n", " res = (t - Phi*m_N)\n", " beta = (N - gamma) / (res*res)\n", " print gamma, alpha, beta" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\t

### γと2αEw(mN)のグラフ

\n", "\t

\n", "\t\tM=4の多項式回帰に対して、図3.16(a)\n", "\t

\n", "\t

\n", "\t\t\n", "\t

\n", "\t

\n", "\t\tのように$\\gammaと2\\alpha E_w(m_N)$のグラフを$\\ln \\alpha$の関係を図化してみます。\n", "\t

\n", "\t

\n", "\t\t形が図3.16とは異なり、曲線の交差する点もα=0.024だとln αの値は負になるのに、図3.16は正の部分で\n", "\t\t交差しています。（何かちがう！）\n", "\t

### 残差と対数エビデンスのグラフ

\n", "\t

\n", "\t\t同様に、M=4の多項式回帰に対して、図3.16(b)\n", "\t

\n", "\t

\n", "\t\t\n", "\t

\n", "\t

\n", "\t\tのように残差と対数エビデンスを$\\ln \\alpha$の関係を図化してみます。\n", "\t

\n", "\t

\n", "\t\tやはり形が、図3.16とは異なります。\n", "\t

### M=4に対する最適な解

\n", "\t

\n", "\t\tM=4の多項式回帰に対する最適なα、βを使って$sin(2 \\pi x)$、観測データ、フィッティング曲線、分散をグラフに表示してみます。\n", "\t

\n", "\t\tM=9の最適解（図1.17）\n", "\t

\n", "\t\t\n", "\t

\n", "\t\tと比べると若干フィッティングが良くありませんが、与えられた点をスムーズに補完しているのがわかります。\n", "\t

\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

## どうしてエビデンス関数が合わないのか

\n", "\t

\n", "\t\tつぎにどうして、図3.16と合わないのか、調べてみます。\n", "\t

\n", "\t

\n", "\t\t図3.16の説明をよく読むと、基底関数が多項式ではなく、ガウス基底関数になっていました。\n", "\t

\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

### 基底関数のチェック

\n", "\t

\n", "\t\tガウス基底関数がどのようなものなのか、図3.1を再現しながら、見てみましょう。\n", "\t

\n", "\t

\n", "\t\t以下に、\n", "\t\t

\n", "\t\t\t
• 多項式の基底関数
• \n", "\t\t\t
• ガウス基底関数
• \n", "$$\n", "\t\t\t\\phi_j(x) = exp \\left \\{ - \\frac{(x - \\mu_j)^2}{2 s^2} \\right \\}\n", "$$\n", "\t\t\t
• シグモイド基底関数
• \n", "$$\n", "\t\t\t\\phi_j(x) = \\sigma \\left( \\frac{x - \\mu_j}{s} \\right)\n", "$$\n", "$$\n", "\t\t\t\\sigma(a) = \\frac{1}{1 + exp(-a)}\n", "$$\n", "\t\t
\n", "\t\tを表示します。\n", "\t

\n", "\t

\n", "\t\tガウス基底関数、シグモイド基底関数では、μの値として-1から1（問題によって0から1の場合もある）に\n", "\t\tとり、それを基底関数の数で等分した値をとるみたいです。\n", "\t\t線形回帰モデルとかを参考にしました。\n", "\t

### ガウス基底関数を使ったγと2αEw(mN)のグラフ

\n", "\t

\n", "\t\tM=9個のガウス基底関数を使ったγと2αEw(mN)のグラフを以下に示します。\n", "\t

\n", "\t

\n", "\t\t関数の形状と交差点の位置が図3.16と合います。\n", "\t