宣言どおり，
FIT申し込み締め切り近し - なぜか数学者にはワイン好きが多い
Mathematicaによる粒子フィルタの実装結果を公開します．
まず最初に，正解というか，いわゆる「真の値」はこれです．これをノイズが入った観測値から推測するわけです．
そして，下記のプログラムを実行すると，
Timing[plist={};rlist={};clist={};tlist={}; // あとでプロットするための変数
Block[{np=100, nd=100}, // 粒子の数と，ループ回数
pv[m_,sig_]:=RandomReal[NormalDistribution[m,sig]]; // 正規分布による乱数
pf[x_]:=x+ pv[0,0.1]; // システムモデル
ph[x_]:=x+ pv[0,0.1]; // 観測モデル
trueval[k_]:=Which[k<25,1,k<50,0,k<75,-1,k<=100,1]; // 観測できない真の値
alpha[y_,p_,sig_]:=1/(Sqrt[2 Pi ]sig)E^(-((y-p)^2/(2 sig^2))); // 粒子の尤度
newfs=Table[pf[trueval[0]],{np}]; // 再サンプリングの初期値
Do[fs=Table[pf[newfsi],{i,1,np}]; // 一期前の粒子からの一期先予測の計算
y=ph[trueval[j]]; // 観測値の取得
alphalist=Table[alpha[y,fsi,0.1],{i,1,np}]; // 粒子の尤度の計算
totalalpha=Total[alphalist]; // 正規化定数
newfs=RandomChoice[alphalist/totalalpha->fs,np]; // 再サンプリング
plist=Append[plist,Mean[newfs]]; // 再サンプリング結果の平均の記録
rlist=Append[rlist,Mean[fs]]; // 一期先予測結果の平均の記録
clist=Append[clist,y]; // 観測値の記録
tlist=Append[tlist,trueval[j]] // 真の値の記録
,{j,1,nd}]]]
次のような結果が得られます．青がノイズの入った観測値，赤がノイズの入ったシステムモデルによる予測値，黒が再サンプリングをした最終的な予測値です．
システムノイズがガウス分布で分散が小さいので，余り観測値を信用せずに自分のシステムモデルの方を信用するために，観測値を追うのが遅れています．
じゃあ，分散を大きくするとどうなるかと言いますと，
pf[x_] := x + pv[0, 1] // システムモデル
観測値を信用してくっついていますが，そのぶん，ガタガタになってます．
さて，ではコーシー分布を用いた場合ですが，Mathematicaのプログラムを次のようにちょっとだけ改良します．
ps[m_, sig_] := RandomReal[CauchyDistribution[m, sig]]; // コーシー分布によるノイズ
pv[m_, sig_] := RandomReal[NormalDistribution[m, sig]]; // ガウス分布によるノイズ
pf[x_] := x + ps[0, 0.02]; // コーシー分布によるシステムモデルのノイズ
すると，次のようになります．
赤い線は，コーシー分布により生成された粒子で，ときたま激しく観測値と離れた値をとります．しかし，再サンプリングにより補正されて，黒い線で表されるように，結果的に非常に「真の値」に近い予測をしています．
いかがでしょう．
こんな短いプログラムで，「21世紀の統計科学III」p. 328のグラフと，ほぼ同じ図が得られました．
なお，このプログラムの検証には，次の本も参考にさせて頂きました．お礼を申し上げます．
