Julia によるかんたんなパーティクルフィルタの例

ローカルレベルモデルです。

（システムモデル）
（観測モデル）

パーティクルフィルタで固定パラメータの とか を推定する方法は色々提案されてたりするんだけど、シミュレーテッドアニーリング法に丸投げです。

「こう書くともっと速くなるよ」みたいなご指摘をお待ちしております。

using Distributions
using Random
using StatsFuns
using Optim
using RCall
using Plots

function locallevel(y,N_particles,mu,Sigma,Tau)
  rng = Random.default_rng()
  tmax = length(y)
  xx = zeros(tmax+1,N_particles)
  xx[1,1:N_particles] = rand(rng, Normal(mu,Tau),N_particles)
  ll = 0
  for t in 2:(tmax+1)
    x = xx[t-1,1:N_particles] + rand(rng, Normal(0,Tau),N_particles)
    w = [logpdf(Normal(x[i],Sigma),y[t-1]) for i in 1:N_particles]
    ll += logsumexp(w) - log(N_particles)
    x = wsample(rng,x,softmax(w),N_particles)
    xx[t,:] = x
  end
    return xx, ll
end

Nile = rcopy(R"Nile")

function loglik(par,y)
  state, lp = locallevel(y,2500,par[1],exp(par[2]),exp(par[3]))
  return -lp
end

Random.seed!(1)
opt = optimize(par -> loglik(par,Nile),[Nile[1],2.0,2.0], SimulatedAnnealing(),Optim.Options(show_trace=true, iterations = 1000))
println(opt)
rawpar = Optim.minimizer(opt)
state, lp = locallevel(Nile,2500,rawpar[1],exp(rawpar[2]),exp(rawpar[3]))
sbar = mean(state[2:101,:],dims=2)


plot(Nile,legend=false)
plot!(sbar,legend=false)
png("./Desktop/Nile.png")