Rで状態空間モデル（MCMC実験編）

先日のエントリーではRでマルコフ連鎖モンテカルロ（MCMC）を実行するための環境構築を行いました。

Rで状態空間モデル（環境構築編）

統計分析、特に時系列分析についていろいろな手法を調べてると、状態空間モデルに辿り着くことが多いと思います。状態空間モデルって何？という方はこちらの素晴らしいサイトをご参照ください。（調べている間に何度も登場したはずです。）logi...

今回は架空のデータを生成してMCMCを実行し、モデルのパラメータをちゃんと推定できるかを実験してみたいと思います。

毎度のことながら馬場先生のサイトを参考にさせていただきました。

logics-of-blue.com

想定するモデルと架空データの生成
パラメータの推定
結果の確認
終わりに

想定するモデルと架空データの生成

まずは、下記の式に従った完全な架空のデータを生成します。

$\\y_t=\beta_{0,t}+\beta_{1,t}x_{1,t}+\beta_{2,t}x_{2,t}+\epsilon_t\\ \beta_{0,t}=\beta_{0,t-1}+\epsilon_{0,t}\\ \beta_{1,t}=\beta_{1,t-1}+\epsilon_{1,t}\\ \beta_{2,t}=\beta_{2,t-1}+\epsilon_{2,t}\\$

つまり２つの時系列データ $x_1$ 、 $x_2$ による線形式から生成される時系列データ $y$ を対象とします。

係数の $\beta$ はランダムウォークする時変係数なので、これはいわゆる動的線形モデル（DLM）ですね。

これらは、観測できないデータ、つまり状態変数です。

さらに今回は $x_1$ 、 $x_2$ も架空データとしてランダムウォークする値を生成しました。

コードはこちら

# パラメタの設定
# サンプルサイズ
N <- 200
# 観測誤差の標準偏差
observationErrorSD <- 20
# 切片の過程誤差の標準偏差
processErrorSD <- 10
# 係数1の過程誤差の標準偏差
coefErrorSD1 <- 0.5
# 係数2の過程誤差の標準偏差
coefErrorSD2 <- 0.2
# 各種シミュレーションデータの生成
# 説明変数をシミュレーションで作る
set.seed(1)
explanatory1 <- rnorm(n=N, mean=10, sd=10)
set.seed(2)
explanatory2 <- rnorm(n=N, mean=20, sd=5)
# slopeのシミュレーション
set.seed(12)
slope1 <- cumsum(rnorm(n=N, mean=0, sd=coefErrorSD1)) + 10
set.seed(13)
slope2 <- cumsum(rnorm(n=N, mean=0, sd=coefErrorSD2)) + 5
plot(slope1, main="時間によるslope1の変化", xlab="day")
plot(slope2, main="時間によるslope2の変化", xlab="day")
# interceptのシミュレーション
set.seed(3)
intercept <- cumsum(rnorm(n=N, mean=0, sd=processErrorSD))
plot(intercept, main="時間によるinterceptの変化", xlab="day")
# responseのシミュレーション
set.seed(4)
response <- intercept + explanatory1*slope1 + explanatory2*slope2 + rnorm(n=N, mean=0, sd=observationErrorSD)
plot(response, main="responseのシミュレーション結果", xlab="day")

ランダムウォークするパラメータのグラフはこちら。

生成された $y$ はこちらとなります。なんだかよくわかりません。

パラメータの推定

さて、これら（ $y,x_1,x_2$ ）のデータに対しさっそく冒頭のモデル式を仮定してMCMCによるパラメータ推定を行います。

MCMCを本当に一言で言うと、

推定したいパラメータの事後分布からのサンプルを、事前分布、確率モデル、データからイイ感じに発生させてくれるアルゴリズム

といった感じです。

パラメータのサンプルを発生させてくれるので、計算回数を多くすればその分サンプルが得られ、

その平均値や分散を計算することやヒストグラムを描くこともでき、

パーセンタイルをとることで推定値の信頼区間を求める（これが大きい！）こともできます。

さて、MCMCを実行するためにはRconsole上で操作するだけではなく、別途stanコードを書く必要があります。

以下のスクリプトを.stanという拡張子でRのカレントディレクトリに保存します。

なぜか最終行にはスペースを入れるというのがルールになっているようです。

data {
int n_sample;       // サンプルサイズ
real y[n_sample];   // 観測値
real x1[n_sample];   // 説明変数1
real x2[n_sample];   // 説明変数2
}
parameters {
real beta0_zero;       // 状態の初期値
real beta1_zero;       // 状態の初期値
real beta2_zero;       // 状態の初期値
real beta0[n_sample];  // 状態の推定値
real beta1[n_sample];  // 状態の推定値
real beta2[n_sample];  // 状態の推定値
real s_w0;  // 過程誤差の分散
real s_w1;  // 過程誤差の分散
real s_w2;  // 過程誤差の分散
real s_v;  // 観測誤差の分散
}
model {
// 状態の初期値から最初の時点の状態が得られる
beta0[1] ~ normal(beta0_zero, sqrt(s_w0));
beta1[1] ~ normal(beta1_zero, sqrt(s_w1));
beta2[1] ~ normal(beta2_zero, sqrt(s_w2));
// 状態方程式に従い、状態が遷移する
for(i in 2:n_sample) {
beta0[i] ~ normal(beta0[i-1], sqrt(s_w0));
beta1[i] ~ normal(beta1[i-1], sqrt(s_w1));
beta2[i] ~ normal(beta2[i-1], sqrt(s_w2));
}
// 観測方程式に従い、観測値が得られる
for(i in 1:n_sample) {
y[i] ~ normal(beta0[i] + beta1[i]*x1[i] + beta2[i]*x2[i], sqrt(s_v));
}
}

今回は架空のデータを自分で作ったのでモデル式はそれと同じものを書けばいいだけです。簡単ですね。

Rに移って以下のコードを実行します。

#パッケージを呼び出す。
library(rstan)
library(ggplot2)
#計算を早くする
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
# データの準備としてリスト形式で保存
data_stan <- list(
y = response,
x1 = explanatory1,
x2 = explanatory2,
n_sample = N
)
# パラメータの推定
fit_stan_count <- stan(
file = "dlm_2factor.stan",
data = data_stan,
iter = 8000,
thin = 10,
chains = 4,
seed = 1
)

こちら、けっこう時間がかかります。

シンプルなモデルですが私のPC（低スペック）で１時間程度かかりました。

推定のコードを少し説明すると、

iterは何回乱数を発生させるか、thinは間引きの間隔、chainsはiterを何回繰り返すか、です。

最初のほうの不安定な乱数を捨てる個数をwarmupで指定するのですが、今回は指定していません。

ただ、ドキュメントを見るとwarmup = floor(iter/2)と記載があるので4000になっていると思います。

したがって今回の場合は(8000-4000)/10*4=1600がMCMCから得られるサンプルの数になります。

有効なサンプルはこれよりも少なくなるのですが、その数は結果のn_effで確認できます。

seedはデータ生成のコードで登場したのと同じ意味で、特定の乱数列を指定しています。

結果の確認

結果を見てみましょう。

fit_stan_count
>Inference for Stan model: dlm_2factor.
>4 chains, each with iter=8000; warmup=4000; thin=10;
>post-warmup draws per chain=400, total post-warmup draws=1600.
>
>              mean se_mean    sd     2.5%      25%     50%     75%   97.5%   n_eff Rhat
>beta0_zero  -15.37    0.70 25.95   -65.61   -33.02  -15.87    2.09   37.49  1389 1.00
>beta1_zero    8.08    0.03  1.16     5.83     7.31    8.05    8.84   10.50  1559 1.00
>beta2_zero    6.16    0.03  1.16     3.91     5.39    6.17    6.91    8.53  1406 1.00
>beta0[1]    -15.60    0.65 23.94   -60.93   -32.09  -16.16    0.82   32.61  1351 1.00
>beta0[2]    -15.18    0.65 23.88   -60.79   -31.41  -15.48    0.34   34.14  1337 1.00
>beta0[3]    -13.01    0.69 23.95   -56.37   -29.22  -13.23    2.21   35.15  1203 1.00
（省略）
>beta2[198]    4.41    0.05  1.30     1.96     3.53    4.34    5.19    7.21   735 1.00
>beta2[199]    4.38    0.05  1.35     1.89     3.50    4.32    5.18    7.34   759 1.00
>beta2[200]    4.42    0.05  1.38     1.85     3.53    4.34    5.27    7.54   743 1.00
>s_w0         74.03    1.86 39.13    12.75    45.95   68.59   94.24  171.06   441 1.01
>s_w1          0.20    0.01  0.11     0.04     0.12    0.18    0.25    0.47   455 1.00
>s_w2          0.11    0.00  0.07     0.01     0.05    0.09    0.14    0.28   260 1.00
>s_v         406.96    2.39 70.90   279.29   357.59  404.40  450.56  555.47   879 1.00
>lp__       -976.42    4.87 85.53 -1129.32 -1036.27 -981.09 -922.82 -791.99   309 1.01

パラメタが収束しているかを表すRhatですが、目安として1.1以下だとOKだそうです。

肝心の推定値のほうは、下から2～5つ目のmeanを見ましょう。

最初に生成した架空データの標準偏差はこちらでした。

# 観測誤差の標準偏差
observationErrorSD <- 20
# 切片の過程誤差の標準偏差
processErrorSD <- 10
# 係数1の過程誤差の標準偏差
coefErrorSD1 <- 0.5
# 係数2の過程誤差の標準偏差
coefErrorSD2 <- 0.2

stanコードでは分散を推定するようにしたので推定値は上の値の2乗が答えとなりますが、概ね正しく推定できています。

最後に、 $\beta$ の推移を見てみます。

sampling_result <- rstan::extract(fit_stan_count)
model_lambda_smooth <- t(apply(
X = sampling_result$beta0,
MARGIN = 2,
FUN = quantile,
probs=c(0.025, 0.5, 0.975)
))
colnames(model_lambda_smooth) <- c("lwr", "fit", "upr")
# データの整形
stan_df <- cbind(
data.frame(intercept = intercept, time = 1:N),
as.data.frame(model_lambda_smooth)
)
# 図示
ggplot(data = stan_df, aes(x = time, y = intercept)) +
labs(title="サンプリング結果") +
geom_point(alpha = 0.6, size = 0.9) +
geom_line(aes(y = fit), size = 1.2) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3)
sampling_result <- rstan::extract(fit_stan_count)
model_lambda_smooth <- t(apply(
X = sampling_result$beta1,
MARGIN = 2,
FUN = quantile,
probs=c(0.025, 0.5, 0.975)
))
colnames(model_lambda_smooth) <- c("lwr", "fit", "upr")
# データの整形
stan_df <- cbind(
data.frame(beta1 = slope1, time = 1:N),
as.data.frame(model_lambda_smooth)
)
# 図示
ggplot(data = stan_df, aes(x = time, y = beta1)) +
labs(title="サンプリング結果") +
geom_point(alpha = 0.6, size = 0.9) +
geom_line(aes(y = fit), size = 1.2) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3)
sampling_result <- rstan::extract(fit_stan_count)
model_lambda_smooth <- t(apply(
X = sampling_result$beta2,
MARGIN = 2,
FUN = quantile,
probs=c(0.025, 0.5, 0.975)
))
colnames(model_lambda_smooth) <- c("lwr", "fit", "upr")
# データの整形
stan_df <- cbind(
data.frame(beta2 = slope2, time = 1:N),
as.data.frame(model_lambda_smooth)
)
# 図示
ggplot(data = stan_df, aes(x = time, y = beta2)) +
labs(title="サンプリング結果") +
geom_point(alpha = 0.6, size = 0.9) +
geom_line(aes(y = fit), size = 1.2) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3)