データ分析スピリッツ

仕事・生活に役立つ分析ツールを提供。コンテンツ拡充中。

【統計分析】Rで状態空間モデル(MCMC実験編)

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

dataspirits.hatenablog.com

今回は架空のデータを生成して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}\\

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

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

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

さらに今回はx_1x_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")

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

f:id:dataspirits:20200628233539p:plainf:id:dataspirits:20200628233556p:plainf:id:dataspirits:20200628233608p:plain

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

f:id:dataspirits:20200628233740p:plain

パラメータの推定

さて、これら(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<lower=0> s_w0;  // 過程誤差の分散
  real<lower=0> s_w1;  // 過程誤差の分散
  real<lower=0> s_w2;  // 過程誤差の分散
  real<lower=0> 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(低スペック)で1時間程度かかりました。

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

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) 

f:id:dataspirits:20200629015937p:plainf:id:dataspirits:20200629015947p:plainf:id:dataspirits:20200629015956p:plain

うまく推定できているようです。正しい\betaの推移(点)は概ね網掛けの99%信頼区間に入っています。

終わりに

今回は架空のデータを生成して答えがわかっているモデルに対し、MCMCを実行しました。

全体の流れを把握することやコードの書き方、結果がどのように表示されるのかを知ることができたと思います。

実際にモデリングする際は何が正しいパラメータかなんて当然わからないため、答え合わせができません。

データや推定結果からそれらしいものを判断するしかない上、そこへ辿り着くまでの試行錯誤も必要になります。

これらの手法で実データに対しうまくモデリングできた際はまたご紹介させていただきたいと思います。