« [OpenJDK]ソースファイルの文字コードを指定してコンパイルする | トップページ | [R]トレンドのあるモデルの比較(水準成分と傾き成分の推移)(「カルマンフィルタ」(共立出版)pp.76-79) »

2023年10月23日 (月)

[R]ローカルレベルモデルによる体重の長期予測とその95%信頼区間と95%予測区間(「カルマンフィルタ」(共立出版)p.38)

公式ページからダウンロードした入力ファイル「Weight.dat」をカレントディレクトリに置いておくこと。太実線が予測値、細実線が95%信頼区間。破線が95%予測区間。

> dtf <- read.table("Weight.dat", header = FALSE)
> yt <- yt0 <- dtf$V1
> yt <- yt0[1:50]
> n0 <- length(yt0)
> n <- length(yt)
> np <- n0 - n
> t0 <- 1:n0
> t <- 1:n
> tp <- setdiff(t0, t)
> mod <- SSModel(yt ~ SSMtrend(1, Q = NA), H = NA)
> fit <- fitSSM(mod, numeric(2), method = "BFGS")
> kfs <- KFS(fit$model)
> at <- kfs$a[c(-1)]
> ppt <- kfs$P[, , -1] - as.vector(fit$model$Q)
> seta2 <- as.vector(fit$model$Q)
> seps2 <- as.vector(fit$model$H)
> abt <- rep(at[n], np)
> j <- 1:(np)
> ppbt <- rep(ppt[n], np) + j * seta2
> ffbt <- ppbt + seps2
> alp <- 0.05
> abtcon1 <- abt + qnorm(alp / 2) * sqrt(ppbt)
> abtcon2 <- abt + qnorm(1 - alp / 2) * sqrt(ppbt)
> abtpre1 <- abt + qnorm(alp / 2) * sqrt(ffbt)
> abtpre2 <- abt + qnorm(1 - alp / 2) * sqrt(ffbt)
> plot(t0, yt0, type = "n", xlab = "経過日数", ylab = "体重(kg)", ylim = c(83, 87.5))
> lines(t0, yt0, lty = "dotted", col = "gray")
> points(t0, yt0, pch = 21, bg = "white", col = "gray", cex = 1.3)
> lines(tp, abtpre1, lty = "dashed", lwd = 3)
> lines(tp, abtpre2, lty = "dashed", lwd = 3)
> lines(tp, abtcon1, lty = "solid", lwd = 1)
> lines(tp, abtcon2, lty = "solid", lwd = 1)
> lines(tp, abt, lty = "solid", lwd = 2)
> points(t, yt, pch = 21, bg = "white", cex = 1.3)

Kalmanfil_p038

« [OpenJDK]ソースファイルの文字コードを指定してコンパイルする | トップページ | [R]トレンドのあるモデルの比較(水準成分と傾き成分の推移)(「カルマンフィルタ」(共立出版)pp.76-79) »

R(本の計算を再現)」カテゴリの記事

コメント

コメントを書く

コメントは記事投稿者が公開するまで表示されません。

(ウェブ上には掲載しません)

« [OpenJDK]ソースファイルの文字コードを指定してコンパイルする | トップページ | [R]トレンドのあるモデルの比較(水準成分と傾き成分の推移)(「カルマンフィルタ」(共立出版)pp.76-79) »

無料ブログはココログ

■■

■■■