廿TT

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

拡張カルマンフィルタ:ロジスティック方程式

下記のブログで公開されている拡張カルマンフィルタのコードでもうちょっと遊んでみる.

拡張カルマンフィルタによるロトカ・ヴォルテラ方程式へのデータ同化 - 廿TT は当てはまりがいまいちだったので、今度はうまくいく例.

モデルとフィッティング

拡張カルマンフィルタでは, 状態方程式,

{\displaystyle {\boldsymbol {x}}_{t}=f({\boldsymbol {x}}_{t-1})+{\boldsymbol {w}}_{t}}

と観測方程式,

 {\displaystyle {\boldsymbol {z}}_{t}=h({\boldsymbol  {x}}_{{t}})+{\boldsymbol  {v}}_{t}}

を考え, fhヤコビアンを使って状態方程式と観測方程式を線形化する.

ここで {\boldsymbol w}_{t}{\boldsymbol  v}_{t} は正規ノイズである.

ロジスティック方程式

p'=rp(1−p)

を一回差分で近似すると,

r_{t+1}=r_{t}\\
p_{t+1}= r_{t}p_{t}(1-p_{t}) + p_{t}.

これを関数 f(R のコードではGGfunction)として,

p_{t}= p_{t}

を関数 h(R のコードでは FFfunction)とする.

それぞれのヤコビアンは,

 \begin{pmatrix}1 & 0\\
p_t(1-p_t) & r_t-2 p_t r_t+1\end{pmatrix}

と,

 \begin{pmatrix}0  & 1 \end{pmatrix}

である.

以下のカラーテレビ普及率のデータ(第1章第2節3 1 情報通信機器の世帯普及率 : 平成18年版 情報通信白書)にこのモデルの当てはめを行う.

#library(numDeriv)
source("http://www.datall-analyse.nl/EKF.R")
#カラーテレビの普及率
y <- c(0.3, 1.6, 5.4, 13.9, 26.3, 42.3, 61.1, 75.8, 85.9, 90.3,
     93.7, 95.4, 97.7, 97.8, 98.2, 98.5, 98.9, 98.8, 99.2, 99.1,
     98.9, 98.7, 99.0, 99.3, 99.4, 99.3, 99.0, 99.1, 99.0, 98.9,
     99.1, 99.2, 99.2, 98.9, 99.0, 99.2, 99.3, 99.4, 99.0, 99.3, 99.4)/100
#尤度関数
llikss <- function (par, data) {
  mod <- list(
    m0=c(0, 0),#初期状態の期待値
    C0=diag(c(100, 100)), #初期状態の分散
    V=par, #Measurement noise (to be estimated)
    W=diag(c(0,0)) #Process noise
  )
  GGfunction <- function (x, k){
    r <- x[1]
    p <- x[2]
    c(r,
      r*p*(1-p)+p)
  }
  GGjacobian <- function(x,k){
    r <- x[1]
    p <- x[2]
    c(1,0,
      p-p^2,r-2*p*r+1)
  }
  FFfunction <- function (x, k){
    r <- x[1]
    p <- x[2]
    c(p)
  }
  FFjacobian <- function (x, k){
    r <- x[1]
    p <- x[2]
    c(0,1)
  } 
  dlmExtFilter(y=data, mod=mod,
               FFfunction=FFfunction, GGfunction=GGfunction,
               FFjacobian=FFjacobian, GGjacobian=GGjacobian,
               simplify=TRUE, logLik=TRUE)$logLik
}
#最尤推定
opt1 <- optim(par=1, fn=llikss, lower = 0, upper = 10, data=y, method="Brent")
#フィルタリング
mod <- list(
  m0=c(0, 0),
  C0=diag(c(100, 100)),
  V=opt1$par, 
  W=diag(c(0,0)) 
)
GGfunction <- function (x, k){
  r <- x[1]
  p <- x[2]
  c(r,
    r*p*(1-p)+p)
}
GGjacobian <- function(x,k){
  r <- x[1]
  p <- x[2]
  c(1,0,
    p-p^2,r-2*p*r+1)
}

FFfunction <- function (x, k){
  r <- x[1]
  p <- x[2]
  c(p)
}
FFjacobian <- function (x, k){
  r <- x[1]
  p <- x[2]
  c(0,1)
} 
Filt1 <-dlmExtFilter(y=y, mod=mod,
                     FFfunction=FFfunction, GGfunction=GGfunction,
                     FFjacobian = FFjacobian, GGjacobian = GGjacobian)
Smooth1 <-dlmExtSmooth(Filt1)
plot(y)
lines(Filt1$m[-1,2],col="red")
lines(Smooth1$s[-1,2],col="blue")
legend("bottomright",legend=c("filtered","smoothed"),col=c("red","blue"),lty=1)

f:id:abrahamcow:20160919030050p:plain