はじめに
『ゼロからできるMCMC』の図とサンプルコードをR言語で再現します。本と一緒に読んでください。
この記事は、3章「マルコフ連鎖モンテカルロ法の一般論
」の内容です。この章では、マルコフ連鎖モンテカルロ法の典型例としてランダムウォークについて解説します。
【前の章の内容】
【他の章の内容】
【この章の内容】
この章で利用するパッケージを読み込みます。
# 3章で利用するパッケージ library(tidyverse)
3.1 マルコフ連鎖
マルコフ連鎖とは、$x^{(k+1)}$が1つ前の$x^{(k)}$だけに依存し、それ以前の$x^{(1)}, \cdots, x^{(k-1)}$の影響を受けない関係です。
この節では、マルコフ連鎖の典型例としてランダムウォークを行います。ランダムウォークとは、確率$\frac{1}{2}$で$x^{(k+1)} = x^{(k)} + 1$、または$x^{(k+1)} = x^{(k)} - 1$となるモデルです。
初期値が$x^{(0)} = 0$であれば、次のように実装できます。複数回行ってみましょう。
# 繰り返し回数を指定 K <- 1000 # 確率を指定 p <- 0.5 # シミュレーション回数を指定 n_simu <- 100 # シミュレーション res_df <- tibble() # 初期化 for(s in 1:n_simu) { # ランダムウォーク tmp_df <- tibble( k = 1:K, # 繰り返し番号 x = sample(x = c(1, -1), size = K, replace = TRUE, prob = c(p, 1 - p)), # ステップ幅を生成 y = cumsum(x), # k番目までの合計 simulation = as.factor(s) # シミュレーション番号 ) # 結果を結合 res_df <- rbind(res_df, tmp_df) # 途中経過を表示 #print(paste0("simuration: ", s, ", y = ", tmp_df[["y"]][K])) } # 確認 head(res_df)
## # A tibble: 6 x 4 ## k x y simulation ## <int> <dbl> <dbl> <fct> ## 1 1 -1 -1 1 ## 2 2 1 0 1 ## 3 3 1 1 1 ## 4 4 -1 0 1 ## 5 5 -1 -1 1 ## 6 6 1 0 1
確率pで1を、1 - pで-1をランダムに生成します。cumsum()で、各行において1行目からその行までの和を計算します。
全てのシミュレーション結果をプロットして、推移を確認します。
# 補助線のデータフレームを作成 addline_df <- tibble( k = 1:K, root_k = 2 * sqrt(k) ) # 出力の推移をプロット ggplot() + geom_line(data = res_df, aes(x = k, y = y, color = simulation), alpha = 0.5) + # yの推移 geom_hline(aes(yintercept = 0), linetype = "dashed") + # 期待値 geom_line(data = addline_df, aes(x = k, y = root_k), linetype = "dashed") + # 補助線(上) geom_line(data = addline_df, aes(x = k, y = -root_k), linetype = "dashed") + # 補助線(下) theme(legend.position = "none") + # 凡例 labs(title = "Random Walk", subtitle = paste0("p=", p, ", K=", K, ", simulation:", n_simu))
試行回数$K$が大きくなるにつれて中心(平均値)から離れていくシミュレーション結果もあります。$2 \sqrt{K}$で補助線を引いています。どのくらいの確率で中心から離れるのかが気になった方は、二項分布を参照してください。ランダムウォークとの関係は書いてませんが、二項分布については「二項分布の平均と分散の導出：定義式を利用 - からっぽのしょこ」でまとめています。
3.3 非周期性
3.1節で扱ったランダムウォークは、1回(ステップ)ごとの変化量(ステップ幅)が固定されていました。あるタイミングの値(状態)$x^{(k)}$が更新された後に、$x^{(k)}$と同じ状態に戻るのに必要な回数を考えます。先ほどの例だと、最低でも$x^{(k+1)} = x^{(k)} + 1$、$x^{(k+2)} = x^{(k+1)} - 1 = x^{(k)}$の2ステップ必要です。また他の回数でも可能ですが、偶数回のステップに限られます。元の状態に戻ることができる様々なステップ数の内、最大公約数を周期と呼びます。この場合、周期は2です。
この節では、1ステップごとに$-c$から$c$の連続値をとる一様分布に従う乱数(一様乱数)$\Delta x$によって変化するモデルを考えます。式にすると$x^{(k+1)} = x^{(k)} + \Delta x$です。$\Delta x = 0$のとき、1ステップで元の状態に戻ります(そのままの状態を維持します)。また他のどんなステップ数でも戻ることができ、周期は1です。周期が1であるものを非周期的なマルコフ連鎖と言います。
ではやってみましょう。
# 繰り返し回数を指定 K <- 1000 # 一様分布の幅を指定 c <- 3 # シミュレーション回数を指定 n_simu <- 100 # シミュレーション res_df <- tibble() # 初期化 for(s in 1:n_simu) { # ランダムウォーク tmp_df <- tibble( k = 1:K, # 繰り返し番号 x = runif(n = K, min = -c, max = c), # ステップ幅を生成 y = cumsum(x), # k番目までの合計 simulation = as.factor(s) # シミュレーション番号 ) # 結果を結合 res_df <- rbind(res_df, tmp_df) # 途中経過を表示 #print(paste0("simuration: ", s, ", y = ", tmp_df[["y"]][K])) } # 確認 head(res_df)
## # A tibble: 6 x 4 ## k x y simulation ## <int> <dbl> <dbl> <fct> ## 1 1 -1.29 -1.29 1 ## 2 2 2.69 1.40 1 ## 3 3 0.167 1.57 1 ## 4 4 1.75 3.32 1 ## 5 5 -1.05 2.27 1 ## 6 6 -0.611 1.66 1
sample()を使って1か-1をランダムに生成したのを、runif()で-cからcまでの連続値を生成します。
推移をプロットします。
# 補助線のデータフレームを作成 addline_df <- tibble( k = 1:K, root_k = c * sqrt(k) ) # 出力の推移をプロット ggplot() + geom_line(data = res_df, aes(x = k, y = y, color = simulation), alpha = 0.5) + # yの推移 geom_hline(aes(yintercept = 0), linetype = "dashed") + # 期待値 geom_line(data = addline_df, aes(x = k, y = root_k), linetype = "dashed") + geom_line(data = addline_df, aes(x = k, y = -root_k), linetype = "dashed") + theme(legend.position = "none") + # 凡例 labs(title = "Random Walk", subtitle = paste0("c=", c, ", K=", K, ", simulation:", n_simu))
概ね$c \sqrt{K}$の範囲となるようです(どうやって求めるの？)。
参考文献
おわりに
既約性、詳細釣り合い条件についてはいい例がなかったので保留します。以前読んだ本に載ってたような、でもどの本かも覚えてない。
さーんっ歩進んでっ2歩下がる～
【次の章の内容】
