始める無料で始める

回帰の事前分布を可視化する

前の演習では、体重 \(Y\) を身長 \(X\) で予測するベイズ回帰モデル($Y \sim N(m, s^2)$、平均 $m = a + bX$)の各パラメータ(\(a\), \(b\), $s$)について、10,000 個の samples をシミュレートしました。samples の各に含まれる \(a\), \(b\), \(s\) の値は、事前分布として妥当な回帰シナリオの 1 つを表しています。これらの事前シナリオの広がりを探るために、最初の 12 組の事前パラメータ \(a\), \(b\), \(s\) のそれぞれから、身長と体重のペアを 50 組シミュレートします。

この演習はコースの一部です

RJAGS によるベイズモデリング

コースを見る

演習の手順

  • samples の最初の 12 組の事前パラメータを n = 50複製したデータフレーム prior_simulation を作成します(合計 600 行になります)。
  • prior_simulation の 600 行それぞれについて、以下を実行します。
    • \(N(170, 10^2)\) モデルから height の値をシミュレートします。
    • \(N(a + bX, s^2)\) モデルから weight の値をシミュレートします。ここで \(X\) は身長、\((a, b, s)\) は事前パラメータの組です。
  • これで 12 組のパラメータセットそれぞれについて、heightweight のペアが 50 組シミュレートされました。ggplot() を使って、各パラメータの set ごとに 50 組の散布図を作成しましょう。weight を y 軸に配置することを忘れずに!

実践的なインタラクティブ演習

このサンプルコードを完成させて、この演習に挑戦してみましょう。

# Replicate the first 12 parameter sets 50 times each
prior_scenarios_rep <- bind_rows(replicate(n = ___, expr = samples[1:12, ], simplify = FALSE)) 

# Simulate 50 height & weight data points for each parameter set
prior_simulation <- prior_scenarios_rep %>% 
    mutate(height = rnorm(n = 600, mean = ___, sd = ___)) %>% 
    mutate(weight = rnorm(n = 600, mean = ___, sd = ___))

# Plot the simulated data & regression model for each parameter set
ggplot(prior_simulation, aes(x = ___, y = ___)) + 
    geom_point() + 
    geom_smooth(method = "lm", se = FALSE, size = 0.75) + 
    facet_wrap(~ set)
コードを編集して実行