ベイジアンにとってMCMCは非常に強力な手法です。例えば観測データを元にパラメータの事後分布を推定したいときに使えます。しかし、事後分布の形状によっては十分なサンプルサイズを得るまでに非常に大きな計算量を必要とする欠点があります。具体的には、分布が多峰性を持つ場合と、峰が細長い場合です。PRML(パターン認識と機械学習)では後者に対する解決策としてHamiltonian MCという手法が紹介されていますが、モデルが微分可能でなければならないという限界がありますし、多峰性に対する解決策ではありません。ここではPRMLで触れられていない方法のひとつとしてレプリカ交換法を試してみます。
準備
まず、多峰性を持つ分布としてOld Faithful間欠泉のデータをEMアルゴリズムで混合ガウス分布にしたものを用意しました。1行目で読み込んでいるfaithful.txtは
http://www.math.pku.edu.cn/teachers/jjia/Computational%20Statistics/faithful.txt
ここから持ってきました。スペース区切りのデータを2列の行列として読みます。
M = dlmread('faithful.txt');
G = fitgmdist(M,2);
[X,Y] = meshgrid(linspace(1.5,5.5),linspace(40,100));
Z = pdf(G,[X(:),Y(:)]);
fig = figure(1);
clf(fig);
fig.Position = [100,100,1200,900];
hold on;
scatter(M(:,1),M(:,2),'o');
contour(X,Y,reshape(Z,size(X)));
PRML下巻153ページと同じfigureが得られました。
Metropolis法
準備した混合ガウス分布からMCMCでサンプリングすることで、多峰性に対するMCMCの振る舞いを観察します。まず普通のMetripolis法を試してみました。
コードは以下の通りです。
function sample = Metropolis_MCMC(GMM_model, data, time)
sample = zeros(time,numel(data(1,:)));
sample(1,:) = data(randi(size(data,1)),:);
for t = 2:time
diff = randn(1,2).*[2.0,20.0];
candidate = sample(t-1,:) + diff;
L_now = log(pdf(GMM_model,sample(t-1,:)));
L_next = log(pdf(GMM_model,candidate));
if rand < exp(L_next - L_now)
sample(t,:) = candidate;
else
sample(t,:) = sample(t-1,:);
end
end
end
試すためのコードは以下の通りです。
sample = Metropolis_MCMC(G,M,100000);
sample(1:10000,:) = [];
scatter(sample(:,1),sample(:,2),'rx');
fig = figure(2);
clf(fig);
autocorr(sample(:,1),50);
結果として、正しくサンプリングできているように見えますが、別の峰に移るのが困難であるため自己相関が高くなっています。
レプリカ交換法
別の峰に移るのが困難であるのは、峰の尤度に対して鞍の尤度が低すぎるからです。そこで、混合ガウス分布の対数尤度に一定値を掛ける、言い換えると精度行列に一定値を掛けることで鞍を越えやすいようにした分布を考えます。例えば、最初のコードの4行目を
T_inv = 0.1;
Z = exp(log(pdf(G,[X(:),Y(:)]))*T_inv);
とすると、等高線は
となります。掛けている値を逆温度と呼びます。この分布から通常のMetropolis法でサンプル列を得ると、当然ながら自己相関は低くなります。(散布図は省きます)
ここで、異なる逆温度のサンプル列を同時並行的に得つつ、双方のサンプル列の詳細釣り合い条件を保ったままサンプルを交換することができます。具体的には
レプリカ交換法 - Wikipediaなどにありますが、交換する確率を
とします。Eは逆温度を考えないときの対数尤度です。逆温度が低い側の尤度が高いならば必ず交換します。
逆温度が低いサンプル列で鞍を越え、逆温度が高いサンプル列と交換することによって自己相関を低くできるわけです。
コードと結果は以下の通りです。
function [sample, rate] = Replica_Exchange_MCMC(
GMM_model, data, time, replica_num, d_temp)
T_inv = 1.0./(d_temp.^(0:replica_num-1));
sample = zeros(time,numel(data(1,:)),replica_num);
for n = 1:replica_num
sample(1,:,n) = data(randi(size(data,1)),:);
end
rate = [0,0];
for t = 2:time
for n = 1:replica_num
diff = randn(1,2).*[2.0,20.0];
candidate = sample(t-1,:,n) + diff;
L_now = log(pdf(GMM_model,sample(t-1,:,n)));
L_next = log(pdf(GMM_model,candidate));
if rand < exp((L_next - L_now)*T_inv(n))
sample(t,:,n) = candidate;
else
sample(t,:,n) = sample(t-1,:,n);
end
end
for n = 1+rem(t,2):2:replica_num-1
L_lowtemp = log(pdf(GMM_model,sample(t,:,n)));
L_hightemp = log(pdf(GMM_model,sample(t,:,n+1)));
if rand < exp((L_hightemp - L_lowtemp)*(T_inv(n)-T_inv(n+1)))
tmp = sample(t,:,n);
sample(t,:,n) = sample(t,:,n+1);
sample(t,:,n+1) = tmp;
rate(1) = rate(1) + 1;
else
rate(2) = rate(2) + 1;
end
end
end
sample = sample(:,:,1);
rate = rate(1) / sum(rate);
end
sample = Replica_Exchange_MCMC(G,M,100000,6,1.6);
sample(1:10000,:) = [];
scatter(sample(:,1),sample(:,2),'rx');
fig = figure(2);
clf(fig);
autocorr(sample(:,1),50);
