読者です 読者をやめる 読者になる 読者になる

廿TT

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

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

R 微分方程式

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

拡張カルマンフィルタによるロトカ・ヴォルテラ方程式へのデータ同化 - 廿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

拡張カルマンフィルタによるロトカ・ヴォルテラ方程式へのデータ同化

R 微分方程式

※ぼくは「データ同化」のなんたるかについてはまるでわかっていない.

下記のブログで解説されている, 拡張カルマンフィルタを試してみる.

ロトカ・ヴォルテラ方程式については下記エントリを参照されたし.

モデルとフィッティング

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

{\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} は正規ノイズである.

ロトカ・ヴォルテラの方程式,

\frac {dx}{dt}=\alpha x-\beta xy \\
\frac {dy}{dt}=\delta xy-\gamma y

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

x_{t+1}=\alpha x_{t}-\beta x_{t}y_{t} + x_{t} \\
y_{t+1}=\delta x_{t}y_{t}-\gamma y_{t} + y_{t}.

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

x_{t}= x_{t} \\
y_{t}=y_{t}.

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

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

 \begin{pmatrix}1 + \alpha - \beta x & -\beta y \\
\delta y & 1 + \delta x - \gamma \end{pmatrix}

と,

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

である.

ロトカ・ヴォルテラの方程式のパラメータ推定 - 廿TT と同じく, 以下のヤマネコとウサギのデータセットにこのモデルの当てはめを行う.

Hare=c(30, 47.2, 70.2, 77.4, 36.3, 20.6, 18.1, 21.4, 22, 25.4, 27.1, 40.3, 57, 76.6,
       52.3, 19.5, 11.2, 7.6, 14.6, 16.2, 24.7)
Lynx=c(4, 6.1, 9.8, 35.2, 59.4, 41.7, 19, 13, 8.3, 9.1, 7.4, 8, 12.3, 19.5,
       45.7, 51.1, 29.7, 15.8, 9.7, 10.1, 8.6)

まずは最尤法で未知パラメータ(α, β, γ, δ と観測ノイズの分散パラメータ.
今回, システムノイズは 0 とした)の推定を行う.

optim の初期値は, ロトカ・ヴォルテラの方程式のパラメータ推定 - 廿TT でやったのと同じ方法で決めた.

library(cowplot)
library(dlm)
source("http://www.datall-analyse.nl/EKF.R")
data1 <- cbind(Hare,Lynx)
#最尤推定
llikss <- function (par, data) {
  mod <- list(
    m0=c(0, 0),#初期状態の期待値
    C0=diag(c(100, 100)), #初期状態のぶんさん
    V=diag(exp(par[1:2])), #Measurement noise (to be estimated)
    W=diag(c(0,0)) #Process noise
  )
  alpha <-exp(par[3])
  beta <- exp(par[4])
  delta <- exp(par[5])
  gamma <- exp(par[6])
  #WARNING: always use arguments x and k when specifying the GGfunction
  GGfunction <- function (x, k){
    prey <- x[1]
    pred <- x[2]
    c(prey + (alpha*prey - beta*prey*pred),
      pred + (delta*prey*pred - gamma*pred))
  }
  #WARNING: always use arguments x and k when specifying the GGjacobian
  GGjacobian <- function (x, k){
    prey <- x[1]
    pred <- x[2]
    c(1 + alpha - beta*pred, -beta*prey,
      delta*pred, 1 + delta*prey - gamma)
  }
  #WARNING: always use arguments x and k when specifying the FFfunction
  FFfunction <- function (x, k){
    prey <- x[1]
    pred <- x[2]
    c(prey, pred)
  }
  #WARNING: always use arguments x and k when specifying the FFjacobian
  FFjacobian <- function (x, k){
    prey <- x[1]
    pred <- x[2]
    c(1, 0,
      0, 1)}
  dlmExtFilter(y=data, mod=mod,
               FFfunction=FFfunction, GGfunction=GGfunction,
               GGjacobian=GGjacobian, FFjacobian=FFjacobian,
               simplify=TRUE, logLik=TRUE)$logLik
}
dydty=(1/Lynx[2:(n-1)])*((Lynx[3:n]-Lynx[1:(n-2)])/2)
dxdtx=(1/Hare[2:(n-1)])*((Hare[3:n]-Hare[1:(n-2)])/2)

fit1 <-lm(dydty~Hare[2:(n-1)])
fit2 <-lm(dxdtx~Lynx[2:(n-1)])

alpha = unname(-fit1$coefficients[1])
beta = unname(fit1$coefficients[2])
gamma = unname(fit2$coefficients[1])
delta = unname(-fit2$coefficients[2])
parini <- c(var(Hare),var(Lynx),alpha,beta,delta,gamma)
opt1 <- optim(log(parini), llikss, data=data1, method="SANN",control = list(temp=10))
#フィルタリング
pars <- exp(opt1$par)
mod <- list(
  m0=c(0, 0),
  C0=diag(c(100, 100)),
  V=diag(pars[1:2]), 
  W=diag(c(0,0)) 
)
GGfunction <- function (x, k){
  prey <- x[1]
  pred <- x[2]
  c(prey + (pars[3]*prey - pars[4]*prey*pred),
    pred + (pars[5]*prey*pred - pars[6]*pred))
}
GGjacobian <- function (x, k){
  prey <- x[1]
  pred <- x[2]
  c(1 + pars[3] - pars[4]*pred, -pars[4]*prey,
    pars[5]*pred, 1 + pars[5]*prey - pars[6])
}
FFfunction <- function (x, k){
  prey <- x[1]
  pred <- x[2]
  c(prey, pred)
}
FFjacobian <- function (x, k){
  prey <- x[1]
  pred <- x[2]
  c(1, 0,
    0, 1)
}
Filt1 <-dlmExtFilter(y=data1, mod=mod,
                     FFfunction=FFfunction, GGfunction=GGfunction,
                     GGjacobian=GGjacobian, FFjacobian=FFjacobian)
n <-nrow(data1)
p1 <-ggplot()+
  geom_point(aes(x=1:n,y=Hare),colour="red")+
  geom_line(aes(x=1:n,y=Filt1$m[-1,1]),colour="red")+
  geom_point(aes(x=1:n,y=Lynx),colour="blue")+
  geom_line(aes(x=1:n,y=Filt1$m[-1,2]),colour="blue")+
  xlab("")+ylab("")
print(p1)

#スムージング
Smooth1 <-dlmExtSmooth(Filt1)
p1 <-ggplot()+
  geom_point(aes(x=1:n,y=Hare),colour="red")+
  geom_line(aes(x=1:n,y=Smooth1$s[-1,1]),colour="red")+
  geom_point(aes(x=1:n,y=Lynx),colour="blue")+
  geom_line(aes(x=1:n,y=Smooth1$s[-1,2]),colour="blue")+
  xlab("")+ylab("")
print(p1)
#####
library(deSolve)
LVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dx = alpha*x-beta*x*y
    dy = delta*x*y-gamma*y
    return(list(c(dx,dy)))
  })
}
ini=Filt1$m[2,]
names(ini) <- c("x","y")
tim <- seq(1,n,by=1)
out1   <- ode(ini,tim, LVmod, c(alpha=pars[3],beta=pars[4],delta=pars[5],gamma=pars[6]))
p1 <-ggplot()+
  geom_point(aes(x=1:n,y=Hare),colour="red")+
  geom_line(aes(x=1:n,y=out1[,"x"]),colour="red")+
  geom_point(aes(x=1:n,y=Lynx),colour="blue")+
  geom_line(aes(x=1:n,y=out1[,"y"]),colour="blue")+
  xlab("")+ylab("")
print(p1)

フィルタリングの結果をプロットした図がこちら. 点が観測値, 折れ線が状態の推定値である.

f:id:abrahamcow:20160918165121p:plain

スムージングの結果がこちら.

f:id:abrahamcow:20160918165254p:plain

最後に, 最尤法で求めたパラメータを使って, 微分方程式を数値的に解いた結果がこちら.

f:id:abrahamcow:20160918165325p:plain

当てはまりはいまいちだなー.

拡張カルマンフィルタ:ロジスティック方程式 - 廿TT に続く.