生活リズム...

ハルヒの映画を見に行く前に一眠り…と思ったら眠れないでござる。
生活リズム崩壊してもうダメダメね…。

それはそうとPRMLようやく3章まで終わった。と言っても3章の後半何やってるかは(多分)把握したが、ちゃんと手を動かしてないので理解が不十分な気分。あとで読み直そう。

眠れないのでmatplotlibの練習がてらに3.3.1の"線形モデルy(x,w)=w0+w1xに対する逐次ベイズ学習の例"を実際に自分でもやってみた。グラフにラベルとか付けようとしたあたりでやる気がなくなったので中途半端感がすさまじい…。
一応上から観測データが0,1,2,10個の時のグラフ。左端が尤度関数、中央がwの事後分布、右端がwの事後分布から適当に取ってきて回帰関数引いたやつ
…ってまんま本の通りだけど。

#!/usr/bin/python
# -*- coding: utf-8 -*-

from scipy import *
from scipy import linalg
from scipy import stats
import matplotlib.pyplot as plt

a = 2.0
b = 25.0

a0 = -0.3
a1 = 0.5

N = 10

# 線を引く本数
L= 5

def norm2D(x, y, u, s):
    sinv = linalg.inv(s)
    d = len(s)
    x_ = x-u[0]
    y_ = y-u[1]
    return exp(-((x_*sinv[0,0]+y_*sinv[1,0])*x_+(x_*sinv[0,1]+y_*sinv[1,1])*y_)/2)/(2*pi)**(d/2)*sqrt(linalg.det(sinv))

def wahu(x, t, m0, s0, b):
    n = len(x)
    A = zeros((n, 2), dtype=float)
    for i in range(n):
        for j in range(2):
            A[i, j] = x[i]**j

    #    sinv = a*identity(2, dtype=float) +
    s0inv = linalg.inv(s0)
    sninv = s0inv + b*dot(A.T, A)
    sn = linalg.inv(sninv)
    mn = dot(sn, dot(s0inv, m0) + b * dot(A.T, t))
    return mn, sn

if __name__ == '__main__':

    x = random.rand(N)*2-1
    t = a0+a1*x+random.normal(0, 0.2)

    x0 = arange(-1, 1, 0.01)
    y0 = arange(-1, 1, 0.01)

    X, Y = meshgrid(x0, y0)

    fig = plt.figure()
    l = [0, 1, 2, 10]

    s = identity(2, dtype=float)/a
    m = [0, 0]
    j = 1

    for i in range(N+1):
        if i > 0:
            m, s = wahu([x[i-1]], [t[i-1]], m, s, b)
        if i in l:
            if i > 0:
                axel = fig.add_subplot(4,3,j)
                Z = stats.norm.pdf(t[i-1], X+Y*x[i-1], 0.2)
                axel.imshow(Z, origin='lower', extent=[-1.0, 1.0, -1.0, 1.0])

            axem = fig.add_subplot(4,3,j+1)
            Z = norm2D(X, Y, m, s)
            axem.imshow(Z, origin='lower', extent=[-1.0, 1.0, -1.0, 1.0])
            axer = fig.add_subplot(4,3,j+2)
            if i > 0:
                axer.scatter(x[:i], t[:i])
            for k in range(L):
                w = random.multivariate_normal(m, s)
                axer.plot(x0, w[0]+w[1]*x0)
                axer.axis([-1.0, 1.0, -1.0, 1.0])
            j += 3
    plt.show()

眠たかった…とか自分に言い訳しないでもっとマシなコード書けるようにならにゃいと…