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

廿TT

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

岩波データサイエンス Vol.1 の年輪の例題を dlm でやる

DLM

ノイズが正規分布し、変数間の関係が線形の状態空間モデルは動的線形モデル(Dynamic Linear Model; DLM)と呼ばれる。

動的線形モデルは以下の式で表現できる。

y_t =F_t \theta_t +v_t, v_t \sim N(0,Vt)
\theta _t=G_t \theta_{t−1}+w_t, w_t\sim N(0,Wt)

推定問題

現時刻を k とするとき、推定問題はつぎの三つに分類できる。

  • 予測(prediction):時刻(k-n)までのデータに基づいて現在の値 y_k を推定する。
  • フィルタリング(filtering):時刻 k までのデータに基づいて現在の値 y_k を推定する。
  • 平滑(smoothing):時刻(k+n)までのデータに基づいて現在の値 y_k を推定する。

ローカルレベルモデル

ローカルレベルモデルは最初の式で

\theta_t=\alpha_tF_t=1G_t=1v_t = \varepsilon _tV_t=\sigma^{2}_{\varepsilon}w_t=\eta_tW_t=\sigma^2_{\eta}

の場合に相当する。

予測、フィルタリング、平滑を実行する R のコードは以下のようになる。

X <- 1961:1990
Y <- c(4.71, 7.70, 7.97, 8.35, 5.70,
       7.33, 3.10, 4.98, 3.75, 3.35,
       1.84, 3.28, 2.77, 2.72, 2.54,
       3.23, 2.45, 1.90, 2.56, 2.12,
       1.78, 3.18, 2.64, 1.86, 1.69,
       0.81, 1.02, 1.40, 1.31, 1.57)

N <- length(Y)
n <- N-6
Y2 <- Y[1:n]
X2 <- X[1:n]
X3 <- X[(n+1):N]

build1 <- function(theta){
  dlm(FF=matrix(1),
      V=matrix(exp(theta[1])),
      GG=matrix(1),
      W=matrix(exp(theta[2])),
      m0 = c(0),
      C0 = matrix(1e+7))
}

fit1 <- dlmMLE(Y2, parm=c(log(var(Y2)),1), build1)
model1 <- build1(fit1$par)
Filt1 <- dlmFilter(Y2,model1)
Smooth1 <-dlmSmooth(Filt1)
Fore1 <- dlmForecast(Filt1,6) 
plot(X,Y,type="b",ylim=c(0,9))
lines(X2,dropFirst(Filt1$m),col="red",lwd=2)
lines(X2,dropFirst(Smooth1$s),col="blue",lwd=2)
lines(X3,Fore1$f,col="orange",lwd=2)
legend("topright",c("observed","filtered","smoothed","forecasted"),
       col=c("black","red","blue","orange"),lty=1,lwd=2)

f:id:abrahamcow:20160507205405p:plain

二階のトレンドモデル

二階のトレンドモデルは最初の式で、

\theta_t=\begin{pmatrix}\alpha_t\\ \alpha_{t-1}\end{pmatrix}F_t=\begin{pmatrix}1 & 0\end{pmatrix}G_t=\begin{pmatrix}2 & -1\\ 1&0\end{pmatrix}

v_t = \varepsilon _tV_t=\sigma^{2}_{\varepsilon}w_t=\begin{pmatrix}\eta_t\\0 \end{pmatrix}W_t=\begin{pmatrix}\sigma^2_{\eta}&0\\0&0\end{pmatrix}

の場合に相当する。

予測、フィルタリング、平滑を実行する R のコードは以下のようになる。

build2 <- function(theta){
  dlm(FF=matrix(c(1,0),1,2),
      V=matrix(exp(theta[1])),
      GG=matrix(c(2,-1,
                  1,0),2,2,byrow = TRUE),
      W=matrix(c(exp(theta[2]),0,
                 0,0),2,2,byrow = TRUE),
      m0 = c(0,0),
      C0 = matrix(c(1e+7,0,
                    0,0),2,2,byrow = TRUE))
}

fit2 <- dlmMLE(Y2, parm=c(log(var(Y2)),1), build2)

model2 <- build2(fit2$par)
Filt2 <- dlmFilter(Y2,model2)
Smooth2 <-dlmSmooth(Filt2)
Fore3 <- dlmForecast(Filt2,6) 
plot(X,Y,type="b",ylim=c(0,9))
lines(X2,dropFirst(Filt2$m[,1]),col="red",lwd=2)
lines(X2,dropFirst(Smooth2$s[,1]),col="blue",lwd=2)
lines(X3,Fore3$f[,1],col="orange",lwd=2)
legend("topright",c("observed","filtered","smoothed","forecasted"),
       col=c("black","red","blue","orange"),lty=1,lwd=2)

f:id:abrahamcow:20160507211023p:plain

参考文献

岩波データサイエンス Vol.1

岩波データサイエンス Vol.1

カルマンフィルタの基礎

カルマンフィルタの基礎