今、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ベストアンサー
- 回答日時:
確認ですが、
>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
- 回答日時:
A No.1のKulesです。
警告がなくなったようで安心しました。
>肝心のプログラムが間違っているようで、上手く答えが出てくれません。。。
というのは、エラーが出るということではなく、何か「このような答えが出るはず」というパラメータを入れてもそのような答えが出ないということでしょうか?
残念ながら私の使用しているMatlabにはStatistics Toolboxが入っていないので、
poissrndやpoisspdfが動きません。
したがって、あなたの書いたコードが動くかどうかを確かめる術がありません。
また、MCMCがどのようなものなのかを知らないため、どのように動けば正解なのかもわかりません。
少なくともこのプログラム自体にはエラーはなさそう(poissrnd以降が私の環境で動かないので断言はできませんが)なので、あとはアルゴリズム上の問題かと。
お役に立てず申し訳ないです。
ご回答、ありがとうございました。
どうにか書き変えて上手く動きました!
教えていただいたCtrl A + Ctrl Iも使わせて頂いてます!!!
ありがとうございました。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- 数学 偏微分に関して教えてください。 g(t)=f(tx,ty)とおいたとき、g(t)の3階微分と4階微分 3 2023/06/27 21:04
- Visual Basic(VBA) まとめシートから集計シートへA列のコードが一致したら1行コピーするマクロをネット上で見つけました。こ 1 2022/08/30 14:11
- Visual Basic(VBA) エクセル マクロ(A1:A10)までの中で一番多く出た数字をB10に表示 6 2023/04/25 17:01
- Visual Basic(VBA) エクセル VBA 難しいです 1 2023/02/21 15:39
- Visual Basic(VBA) [Excel VBA] このコードでは行の挿入や行の消去をすると13のエラーが出てしまう。 3 2022/12/09 00:29
- Visual Basic(VBA) Excel VBAの解読について質問があります。 概要は、マクロでチェックボックスにチェックすると日 1 2023/02/10 07:50
- Visual Basic(VBA) VBAの繰り返し処理について教えてください。 3 2022/08/02 13:21
- Visual Basic(VBA) excel2021で実行できないマクロ。どこを直したらいいのか 2 2022/03/28 03:40
- その他(プログラミング・Web制作) pandasでまとめてインデックスを削除するにはどうすればいいですか? たとえば、以下のプログラムで 1 2022/07/31 23:09
- Visual Basic(VBA) A列B列C列 3 2023/04/26 18:11
関連するカテゴリからQ&Aを探す
おすすめ情報
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
python03について。
-
文系のSE志望です。プログラミ...
-
Adobe Premiere Proについて質...
-
Adobe Premiere Proです。 シー...
-
python3について。
-
vba クリップボードクリアにつ...
-
python3について。
-
python3について。
-
REGZAに接続できない(パソコン)
-
pythonでのカーソル移動がずれる
-
そのまま使っただけなのに・・...
-
このURLで広告を出しているのは...
-
Arduinoに関する質問
-
プログラムの起動、利用につい...
-
Pythonのコードエラーについて...
-
pythonにてseleniumを使うも、...
-
Google ColaboでGUI作成
-
HTMLソースが表示のページのも...
-
https://paiza.jp/challenges/5...
-
バッチファイルについて
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
texでグラフを横ならびにしたと...
-
std::vectorのマージ
-
Texで図と表を並べたときのキャ...
-
MCMCを用いたパラメータ推定[Ma...
-
LaTeXの矢印位置変更
-
Excel・Word リサーチ機能を無...
-
特定のPCだけ動作しないVBAマク...
-
エクセルで特定の列が0表示の場...
-
UserForm1.Showでエラーになり...
-
配列数式の解除
-
Excel マクロ VBA プロシー...
-
教えて下さい
-
メッセージボックスのOKボタ...
-
一つのTeratermのマクロで複数...
-
String""から型'Double'への変...
-
VBAでfunctionを利用しようとし...
-
ExcelのVBA。public変数の値が...
-
End Sub が必要です。
-
Excel VBAからAccessマクロを実...
-
TERA TERMを隠す方法
おすすめ情報