ここを参考に
http://d.hatena.ne.jp/teramonagi/20110508/1304821958
真の状態も合わせてプロットしてみた
library(FKF) ## alpha[t+1] = alpha[t] + eta[t], eta[t] ~ N(0, 2) ## y[t] = alpha[t] + eps[t], eps[t] ~ N(0, 1) alpha <- cumsum(2.0 * rnorm(100)) y <- alpha + rnorm(100) dt <- matrix(0) ct <- matrix(0) Zt <- matrix(1) Tt <- matrix(1) a0 <- 0 P0 <- matrix(0.01) # objective function objective <- function(par, ...){ -fkf(HHt = matrix(par[1]), GGt = matrix(par[2]), ...)$logLik } ## Estimate parameters: fit.fkf <- optim(c(HHt = 4, GGt = 1), fn = objective, yt = rbind(y), a0 = a0, P0 = P0, dt = dt, ct = ct, Zt = Zt, Tt = Tt) ## Filter data with estimated parameters: fkf.obj <- fkf(a0, P0, dt, ct, Tt, Zt, HHt = matrix(fit.fkf$par[1]), GGt = matrix(fit.fkf$par[2]), yt = rbind(y)) ## Plot data together with filtered value and real state: plot(y) lines(fkf.obj$att[1, ],col = "blue") lines(alpha,col = "red") legend("bottomright", c("original data", "estimation","real state"), col = c("black", "blue","red"), lty = 1)