機械学習に詳しくなりたいブログ

機械学習や数学について勉強した内容を中心に書きます。100%趣味です。記事は数学的に厳密でなかったり誤りを含んでいるかもしれません。ご指摘頂ければ幸いです。

MAP推定をMCMCで解く

先回までにメトロポリス・ヘイスティングズ法によるMCMCの実装を行ってきました。(参考:マルコフ連鎖モンテカルロ法 - メトロポリス・ヘイスティングズ法(1))これによって様々な分布からサンプリングできるようになったので、今回はMAP推定をMCMCで解いてみたいと思います。MAP推定については、線形回帰をMAP推定で解くに書いています。求めたい係数 \mathbf wの事前分布を導入し、ベイズの定理を使って、訓練データが与えられた時の事後分布を求め、その最大値を解とする、というものでした。MCMCを使えば、解析的に解けない問題でも事後分布さえわかっていればそこからサンプリングしてしまえば良いということですね。

さて、いくつかの前提条件を線形回帰をMAP推定で解くと同様のものとすれば、事後分布は


\displaystyle p(\mathbf w|\mathbf t) \propto \exp\left\{ \frac{ -\| \mathbf t - \mathbf X \mathbf w\|^2 }{ 2 \sigma ^2} -\frac{1}{2}(\mathbf{w} - \mathbf{\mathbf m_0})^T \mathbf S_{0}^{-1}(\mathbf{w} - \mathbf{\mathbf m_0}) \right\} \tag{1}

と表されます。この式からMCMCを使って直接サンプリングし、最大値を求めてみます。訓練データは線形回帰でやってきたものと同じく、正弦波に正規分布の誤差が加わっているものにし、3次の多項式


y=w_{0}+w_{1}x+w_{2}x^{2}+w_{3}x^{3}\tag{2}

で回帰することとします。

早速結果です。 \mathbf wの事後分布は以下のようになりました。 f:id:opabinia2:20180421131809p:plain

そしてこの係数を使った近似結果がこちら。 f:id:opabinia2:20180421131918p:plain

ちゃんと近似できていることが確認できました!、、、と言いたいところなのですが、解析的にMAP推定を解いたときと、微妙に値がズレていて、サンプリング数を増やしてもそのズレは解消されず。訓練データなどは全く同じにしているはずなんですが。どちらかが間違っているんでしょうけど、、、近似結果がそれっぽいのでどちらにミスがあるのか見つけられず。うーん。

ミスがあると書いておいて載せるのは微妙なんですが、一応今回のコードです。提案分布の分散が適切でないと、サンプリングが全然進んでいきません。マルコフ連鎖モンテカルロ法 - メトロポリス・ヘイスティングズ法(1)で既に経験済みのはずなのに、また悩んでしまった。

# MAP推定をMCMC法を使って解く
import numpy as np
import matplotlib.pyplot as plt


# y = w0*x^0+....wM*x^M を、引数xの配列数分求める
def y(w, x, M):
    X = np.empty((M + 1, x.size))
    for i in range(M + 1):
        X[i, :] = x ** i
    return np.dot(w.T, X)


# 提案分布
def q(x):
    mean = np.zeros(M+1)
    cov = 0.5*np.eye(M+1)
    z = np.random.multivariate_normal(mean, cov).T
    return x + z


# サンプルする関数
def p(w):
    return np.exp(-np.dot((t - np.dot(X, w)).T, (t - np.dot(X, w)))/(2*(v**2)) - 0.5*np.dot(np.dot(w.T, np.linalg.inv(s0)), w))


# MCMCでMAP推定値を求める
def MAP_MCMC(x0):
    while True:
        # サンプル候補
        w0 = q(x0)
        # 採択するまでサンプル候補をとる
        if p(w0.reshape(M+1, 1))/p(x0.reshape(M+1, 1)) > np.random.uniform(0, 1):
            return w0


# ランダムシードを固定
np.random.seed(0)
# 多項式の最大べき乗数(x^0+...+x^M)
M = 3
# 訓練データ数
N = 10


# 訓練データの列ベクトル
x = np.linspace(0, 1, N).reshape(N, 1)
# 訓練データtの列ベクトル
t = np.sin(2*np.pi*x.T) + np.random.normal(0, 0.2, N)
t = t.reshape(N, 1)

# 行列Xを作成
X = np.empty((N, M+1))
for i in range(M+1):
    X[:, i] = x.reshape(1, N) ** i


# 事前分布の共分散行列
v0 = 1000
s0 = v0 * np.eye(M+1, M+1)

# 真値に加わっている正規分布ノイズの分散
v = 0.2

# サンプル数
N = 100000

# バーンイン
B = 10000

# MCMCの初期値
x0 = np.zeros(M+1)

samples = np.empty([N, M+1])

# MCMC
for i in range(N):
    x0 = MAP_MCMC(x0)
#  print(i)
    samples[i, :] = x0

# print(np.average(samples[:,0]))
# print(np.average(samples[:,1]))
# print(np.average(samples[:,2]))
# print(np.average(samples[:,3]))
#
# plt.hist(samples[B:,0],bins=100,normed=True)
# plt.show()
# plt.hist(samples[B:,1],bins=100,normed=True)
# plt.show()
# plt.hist(samples[B:,2],bins=100,normed=True)
# plt.show()
# plt.hist(samples[B:,3],bins=100,normed=True)
# plt.show()

w0 = np.average(samples[B:,0])
w1 = np.average(samples[B:,1])
w2 = np.average(samples[B:,2])
w3 = np.average(samples[B:,3])

w = np.array([w0,w1,w2,w3]).reshape(M+1,1)

# 求めた係数wを元に、新たな入力x2に対する予測値yを求める
x2 = np.linspace(0, 1, 100)
y1 = y(w, x2, M)
# 結果の表示
plt.xlim(0.0, 1.0)
plt.ylim(-1.5, 1.5)
plt.scatter(x, t)
plt.plot(x2, y1.T)
plt.show()