春休み

気付いたら春休み。単位は…まぁ多分取れてるでしょう。
サークルの方は休みに入ってから勉強会+読書会ラッシュ。RWH読書会、組合せ最適化読書会、ICPC勉強会、アルゴリズムイントロダクション読書会。他にも僕が参加していないのだとデザパにC#にお絵描きの勉強会といった感じ。


サークル以外の生活では、ついカッとなってPRML読み始めた。機械学習とかパターン認識にちょっと興味が…とは言え学部2回で手を出すなら別の選択肢もあった気がする。*1


何か手を動かしたくなったので1章の曲線フィッティングと3章のベイズ線形回帰の所を自分でも適当にコード書いてみた。後者の方はちゃんと逐次学習してねえじゃねえか! とかコード汚ないとかいうのは僕のせいですはい。。。というかこんなんで良いのかな…


行列の計算とかグラフの描画はscipy+matplotlibを利用したけど、こういうのってRとか使うと楽だったりするのかなぁ。

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

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

N = 9
a = 0.2


def fit(x, t, m):
    n = len(t)
    A = zeros((n, m), dtype=float)
    for i in range(n):
        for j in range(m):
            A[i, j] = x[i]**j
    return dot(pinv(A), t)

def draw(w, label=""):
    x = arange(0, 256.0)/128.0*pi
    y = arange(0, 256.0)*0
    for i in range(len(w)):
        y += w[i]*(x**i)
    plt.plot(x, y, label=label)

if __name__ == '__main__':

    x0 = arange(0, 256.0)/128*pi
    y0 = sin(x0)
    plt.plot(x0, y0, label="sin(x)")

    x = random.rand(N)*2*pi
    t = random.normal(sin(x), a)
    plt.scatter(x, t, label="data")

    draw(fit(x, t, 1), "M=1")
    draw(fit(x, t, 4), "M=4")
    draw(fit(x, t, 10), "M=10")

    plt.xlim(0, 2*pi)
    plt.ylim(-1.5, 1.5)
    plt.legend()
    plt.show()

実行してみると


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

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

a = 0.2
a0 = -0.3
a1 = 0.5
N = 20

# 線を引く本数
L= 10

def wahu(x, t, a, 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) + b*dot(A.T, A)
    s = linalg.inv(sinv)
    m = b*dot(dot(s, A.T), t)
    return m, s

if __name__ == '__main__':

    x = random.rand(20)*2-1
    t = a0+a1*x+random.normal(0, a)
    x0 = arange(-1, 1, 0.1)

    while True:
        try:
            n = input('> ')
        except:
            break
        plt.clf()
        if n > 0:
            plt.scatter(x[:n], t[:n])
        m, s = wahu(x[:n], t[:n], a, (1/a)**2)
        for i in range(L):
            l = random.multivariate_normal(m, s)
            plt.plot(x0, l[0]+l[1]*x0)
        plt.xlabel('x')
        plt.ylabel('y')
        plt.xlim(-1, 1)
        plt.ylim(-1, 1)
        plt.show()

観測したデータが0個の時、wの事後分布からランダムに選んだ回帰関数(*10)

観測したデータが2個の時

観測したデータが4個の時

観測したデータが20個の時

*1:フリーソフトでつくる音声認識システム とか言うのは読んだけど、もう1ステップ別の本入れてもよかった気がするなぁ…