初めてのベイズ統計モデリング実践と失敗集

2021年12月31日 11:33

背景

ベイズ統計モデリングをしたいなと思い、馬場真哉さんの以下の実践Data Scienceシリーズを読みました。

とても素晴らしい本でした。ベイズ統計モデリングを平易に書いてくれていて、実際に始めてみようと思わせてくれました。
読んだ後のファーストステップとして、実データに適用したいと思ったのが始まりです。
ところで、プライベートで友達がPS5を発売初期に買っていて、いろいろなゲームを満喫していたのを見て、私も欲しいなと思っていましたが、この年末に眺めてみると、ネットの実売価格が２倍以上高くなっていました。もちろん諦めました。。。。悔しいです。

データ的にみると、PS5の発売から１年以上経っていてカカクコムの価格推移がかなり貯まっていたので、それを平滑化トレンドモデルに適用しつつ価格遷移の説明変数も組み込めると買い時も分かるのではというのがモチベーションです。

先に結論から言うとモデル化できましたが、説明変数は見つけられず、どりふと成分もないという役に立たない結果になりました。私と同じモデリング初心者の方に、実データに当てはめると最初は失敗するので気にしなくていいよ、と勇気づけられたら嬉しいと思っています。

データ

データはカカクコムのPS5の平均価格の推移を目的変数として、同時期のGoogle トレンドの「PS5」のワードのスコア推移を説明変数にしています。

コードを書くときは、本からのほんの少しの発展として、上記の本の著者さんのページにある収束させるためのコツを参考にしています。

手法

本に則って、データ作成とstanの実行をしていきます。データリストは、

データ数(n = 398)
平均価格の推移: (2020-11-06~2021-12-08)
Google トレンドのスコア推移(2020-11-06~2021-12-08)

を入れています。Google トレンドのデータは週次データになっているので、特定の週は全て同じスコアが入っています。

data_list <- list(
  T = length(ps5_data$mean_price),  // データ数
  y = ps5_ata$mean_price,           // 平均価格
  ps5_trend = ps5_explain_data)     // トレンドのスコア推移

stanのコードは本をベースにしつつ、先ほど添付したサイトのテクニックを利用しています。

data {
  int T;                 // データ取得期間の長さ
  vector[T] y;           // 観測値
  vector[T] ps5_trend;   // PS5のトレンド
}

parameters {
  real<lower=0> s_w;   // 水準成分の変動の大きさを表す標準偏差
  real<lower=0> s_z;   // ドリフト成分の変動の大きさを表す標準偏差
  real<lower=0> s_v;   // 観測誤差の標準偏差
  real<lower=0> s_t1;  // 時変成分1の標準偏差
  
  vector[T] mu_err;    // 水準成分の増減量
  vector[T] delta_err; // ドリフト成分の増減量
  vector[T] b1_err;    // PS5トレンドスコアの説明変数の増減量
}

transformed parameters {
  vector[T] mu;        // 水準+ドリフト成分+説明変数の推定値
  vector[T] delta;     // ドリフト成分の推定値
  vector[T] b1;        // 半導体のトレンドスコア説明変数の推定値
  
  // 1時点目
  mu[1] = mu_err[1];
  delta[1] = delta_err[1];
  b1[1] = b1_err[1];
  
  // 2時点目以降
  for (t in 2:T) {
    mu[t] = mu[t-1] + delta[t-1] + b1[t-1] * ps5_trend[t-1] + s_w * mu_err[t];
    delta[t] = delta[t-1] + s_z * delta_err[t];
    b1[t] = b1[t-1] + s_t1 * b1_err[t];
  }
}

model {
  // 標準偏差1の正規乱数を得る
  mu_err[2:T] ~ normal(0,1);
  delta_err[2:T] ~ normal(0,1);
  b1_err[2:T] ~ normal(0,1);

  // 観測方程式
  y ~ normal(mu, s_v);
}

stanの実行パラメータは以下のようにしています。

local_linear_trend_time_varying_coef_stan_modeling <- stan(
  file = "local-linear-trend_time-varying-coef.stan",
  data = data_list,
  seed = 1,
  control = list(adapt_delta = 0.95, max_treedepth = 15))

controlのパラメータについて、計算終了時の warning の内容に基づいて適宜追加しています。
実行は４、５年前の macbook pro で行っていて、パラメータの推定時間は３時間くらいはかかりました。計算時間を考えると、ベイズモデリングする前にクラシックな手法を使って解決できるなら、その方が良いことがわかります。

結果

結果の確認をしていきます。ほぼ参考図書の手順通りに実行していきます。

print(local_linear_trend_time_varying_coef_stan_modeling, 
      pars = c("s_w", "s_z", "s_v", "s_t1"), 
      probs = c(0.025, 0.5, 0.975))

Inference for Stan model: local-linear-trend_time-varying-coef.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean     sd    2.5%     50%   97.5% n_eff Rhat
s_w  1256.92   26.27 299.23  435.98 1305.51 1707.44   130 1.03
s_z    38.89    3.28  37.66    1.41   26.64  150.97   132 1.03
s_v  2677.63    6.59 141.60 2408.78 2674.68 2949.75   462 1.01
s_t1    2.01    0.14   1.88    0.06    1.38    6.88   179 1.02

Samples were drawn using NUTS(diag_e) at Sun Dec 26 13:57:59 2021.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Rhatを見る限り、そこまで悪くないかなという結果です。引き続き、収束の確認をしていきます。

mcmc_combo(local_linear_trend_time_varying_coef_stan_modeling_v3, 
           pars = c("s_w", "s_z", "s_v", "s_t1"))

目視で見ても変なことにはなっていなさそうです。

続いて、モデル、ドリフト成分、ps5のトレンド係数を図示していきます。図示に利用した plotSSM関数は書籍内で定義されている関数です。github 上にも公開されています。本当に著者の方には感謝です。

mcmc_sample <- rstan::extract(local_linear_trend_time_varying_coef_stan_modeling)
p_all <- plotSSM(
  mcmc_sample = mcmc_sample,
  time_vec = ps5_data$date,
  state_name = "mu",
  graph_title = "Local Linear Trend Modeling With Coef Result",
  y_label = "all")
print(p_all)

p_delta <- plotSSM(
  mcmc_sample = mcmc_sample,
  time_vec = ps5_data$date,
  state_name = "delta",
  graph_title = "Local Linear Trend Modeling With Coef Result",
  y_label = "drift")
print(p_delta)

p_new_ps5 <- plotSSM(
  mcmc_sample = mcmc_sample,
  time_vec = ps5_data$date,
  state_name = "b1",
  graph_title = "Local Linear Trend Modeling With Coef Result",
  y_label = "google_ps5_trend")
print(p_new_ps5)

全体のモデリングはうまくいきました。これを見ると（カカクコムの数値から同然ですが）３月ごろが一番買い時でした。
ドリフト成分について、下降と上昇のトレンドが追いきれていない。。。もう少し精査する必要があるのかなと。
Google トレンドのPS5の係数ですが、全期間において係数の確率分布が正と負の両方に分布しています。ほぼ影響がないように見えます。

まとめ

今回の結果だと、ドリフト成分も説明変数とも関連がなく買い時もよくわからないという結果になりました。
iter のパラメータを増やすと、もう少しうまくいくのかなと思いましたが、そうすると半日以上潰れるので、コードの改善やGPU使って高速化？みたいなところも考えたいなと思います。
そもそも統計モデリングをもっと勉強した方がいい気もしますが。
最後までありがとうございます。

この記事が気に入ったらサポートをしてみませんか？