廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

JuliaCallを試す(パーティクルフィルタ)

今日の川柳

ジュリアで簡単なパーティクルフィルタを書いた。
R によるすごくかんたんなパーティクルフィルタの実装例 - 廿TT

using Distributions
function PF(y,N_particles,Sigma,Tau)
    tmax = length(y)
    xx = zeros(tmax+1,N_particles)
    xx[1,1:N_particles] = y[1]
    for t in 2:(tmax+1)
      x = xx[t-1,1:N_particles] + rand(Normal(0,1),N_particles) .* Tau
      w = exp.((-(y[t-1]-x).^2)./(2*(Sigma.^2))) ./ Sigma
      ind = wsample(1:N_particles,w/sum(w),N_particles)
      x = x[ind];
      xx[t,1:N_particles] = x;
    end
      return xx
end

これを R から呼ぶ。

library(JuliaCall)
julia_source("./Documents/PFtest.jl")
PFout <-julia_call("PF",Nile,10000L,10,10)

うーん、早いんだろうか。

R で書いたのと比べてみる。

particle_filter <- function(y,N) {
  tmax <- length(y)
  xx <- matrix(NA_real_,tmax+1,N)
  xx[1,] <-rep(y[1],N)
  for(t in 1:tmax){
    x <- xx[t,] + rnorm(N,0,10) #一期先予測
    w <- dnorm(y[t],x,10) #尤度
    xx[t+1,] <- sample(x,N,replace=TRUE,prob=w/sum(w)) #リサンプリング
  }
  return(xx)
}
system.time({
out1 <-particle_filter(Nile,1000000L)
})
#  user  system elapsed 
#29.037   2.762  31.932 
system.time({
out2 <-julia_call("PF",Nile,1000000L,10,10)
})
#   user  system elapsed 
# 25.098   2.440  28.112

微差だなあ。もっとうまい書き方があるのかな?

結果は一致しました。

mean1 <-rowMeans(out1[-1,])
mean2 <-rowMeans(out2[-1,])
plot(Nile)
lines(start(Nile)[1]:end(Nile)[1],mean1,col="royalblue",lwd=2)
plot(Nile)
lines(start(Nile)[1]:end(Nile)[1],mean2,col="red",lwd=2)

f:id:abrahamcow:20180721030920p:plain

f:id:abrahamcow:20180721030941p:plain