MCMCを用いたパラメータ推定[Matlab]

質問

今、MCMCを用いてパラメータ推定のためのプログラムを考えています。Matlabを用いてコードを書いている最中なのですが、IfとEndが合わないという忠告が何個も出てしまいます。下記にコードを載せましたので、何が間違っているか教えていただけませんでしょうか。
求めたいパラメータは、ポアッソン分布に従うと仮定し,尤度を使ってデータの当てはまりを見ています。
Matlab初心者なので間違いも多いと思いますが、どうぞよろしくお願いいたします。

logpt =0; %空値
theta =0; %空値
theta0 = ceil(100.*rand()); % i = initial value of MCMC 100
lambda = 5; % true mean value of dispersal
updatewidth = 0.1; %sigma update width of MCMC
nsample = 100; %number of iterations for MCMC nsample 100000
nthin = 1; %thinning number of MCMC kannbiki, sukashi
individual = 100; %individual numbers 1000
xsize = 10; %space size of matrix
ysize = 10;
dat = zeros(xsize, ysize); %space
xinit = 5; %initial position of population
yinit = 5;

%making dammy data
for j = 1:individual
lam(1) = ceil(100.*rand()); % decide a initial number
P = poissrnd(lambda);
tx = xinit;
ty = yinit;
i = 1;
while i <= P
d = ceil(4.*rand());
if d==1 && ty -1 >0
ty = ty-1;
i = i+1;
elseif 1==2 && tx +1 <= xsize
tx = tx+1;
i = i+1;
elseif d==3 && ty+1 <= ysize
ty = ty+1;
i = i+1;
elseif d==4 && tx-1 > 0
tx = tx-1;
i = i+1;
end;
end;
dat(tx,ty) = dat(tx,ty) + (1);
end;

%MCMC start
for n = 1:nsample
for o = 1:nthin
%(ii) simulation with the previous step dispersal value
sp0 = zeros(xsize, ysize);
ind0 = zeros(2,individual);
for m = 1:individual
r = poissrnd(theta0);
tx = xinit;
ty = yinit;
i = 1;
while i <= r
d = ceil(4.*rand());

if d==1 && ty-1 > 0
ty = ty-1;
i = i+1;
elseif d==2 && tx+1 <= xsize
tx = tx+1;
i = i+1;
elseif d==3 && ty+1 <= ysize
ty = ty+1;
i = i+1;
elseif d==4 && tx-1 > 0
tx = tx-1;
i = i+1;
end;

end;
sp0(tx,ty) = sp0(tx,ty) + 1;
ind0(1,m) = tx;
ind0(2,m) = ty;
end;

%calculating likelihood
LOGP = 0;
for q = 1:xsize
for s = 1:ysize
if sp0(q,s) == 0
LOGP = LOGP + log(poisspdf(0.001,dat(q,s)));
else
LOGP = LOGP + log(poisspdf(sp0(q,s),dat(q,s)));
end;
end;
end;

%(iii) making a proposal value
tmp = updatewidth * randn(size(theta0));
%b= updatewidth = 0.1; sigma update width of MCMC

while tmp < 0
tmp = updatewidth * randn(size(theta0));
end;

%(iv) making proposal distances of each individuals
theta1 = poissrnd(tmp);

%(v) simulation with the proposal value
for p = 1:individual
tx = xinit;
ty = yinit;
i = 1;
while i <= theta1
d = ceil(4.*rand());

if d==1 && ty-1>0
ty = ty-1;
i = i+1;
elseif d==2 && tx+1 <= xsize
tx = tx+1;
i = i+1;
elseif d==3 && ty+1 <= ysize
ty = ty+1;
i = i+1;
elseif d==4 && tx-1 > 0
tx = tx-1;
i = i+1;
end;
end;

%calculating likelihood
prev = sp0(ind0(1,p),ind0(2,p));
dat_prev = dat(ind0(1,p),ind0(2,p));
prop = sp0(tx,ty);
dat_prop = dat(tx,ty);

if prop == 0

LOGP0 = log(poisspdf(prev,dat_prev)) + log(poisspdf(0.001,dat_prop));
else
LOGP0 = log(poisspdf(prev,dat_prev)) + log(poisspdf(prop,dat_prop));
end;
if prev == 1
LOGP1 = log(poisspdf(0.001,dat_prev)) + log(poisspdf(prop+1,dat_prop));
else
LOGP1 = log(poisspdf(prev-1,dat_prev)) + log(poisspdf((prop+1),dat_prop));
end;
if rand < exp(LOGP1 - LOGP0)
r = theta1;
LOGP0 = LOGP1;
end;
end;
%(vi) calculating the mean of theta
theta0 = mean(r);
end;
theta = theta0;
LOGPt = LOGP;
end;

回答No.1

Kules · ベストアンサー

確認ですが、
>IfとEndが合わないという忠告が何個も出てしまいます
というのは「警告」であって「エラー」ではないんですよね？
まあ警告ならほっといてもいいんですが…

もし警告だとすれば、それはifとendのインデントがあっていないことによると思われます。

スペースを_（アンダーバー）で表します。
if 1;
______if 0;
________A=5;
____end;
end;
とすると、おそらく2行目に警告が出ると思われます。
どうやら、ifとendがスペース2つ分以上離れていると出るみたいです。

で、解決方法ですが、Matlabのエディタを使っているのであれば
ctrl+A（全て選択）→ctrl+I（スマートインデント）で解決すると思われます。

参考になれば幸いです。

お礼

回答、ありがとうございました。
教えていただいたやり方で試したところ、警告はなくなりました！
Ctrl+AとCtrl+I、便利ですね！

ただ、肝心のプログラムが間違っているようで、上手く答えが出てくれません。。。
このプログラムに対しても、何かアドヴァイス等ございましたら申し訳ありませんがよろしくお願いします。

ご回答、ありがとうございました。

回答No.2

Kules

A No.1のKulesです。

警告がなくなったようで安心しました。

>肝心のプログラムが間違っているようで、上手く答えが出てくれません。。。
というのは、エラーが出るということではなく、何か「このような答えが出るはず」というパラメータを入れてもそのような答えが出ないということでしょうか？

残念ながら私の使用しているMatlabにはStatistics Toolboxが入っていないので、
poissrndやpoisspdfが動きません。
したがって、あなたの書いたコードが動くかどうかを確かめる術がありません。

また、MCMCがどのようなものなのかを知らないため、どのように動けば正解なのかもわかりません。

少なくともこのプログラム自体にはエラーはなさそう（poissrnd以降が私の環境で動かないので断言はできませんが）なので、あとはアルゴリズム上の問題かと。

お役に立てず申し訳ないです。

お礼

ご回答、ありがとうございました。
どうにか書き変えて上手く動きました！
教えていただいたCtrl A + Ctrl Iも使わせて頂いてます！！！
ありがとうございました。