からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

【R】3.1,3:ランダムウォーク【ゼロからMCMCのノート】

はじめに

 『ゼロからできるMCMC』の図とサンプルコードをR言語で再現します。本と一緒に読んでください。

 この記事は、3章「マルコフ連鎖モンテカルロ法の一般論 」の内容です。この章では、マルコフ連鎖モンテカルロ法の典型例としてランダムウォークについて解説します。

【前の章の内容】

www.anarchive-beta.com

【他の章の内容】

www.anarchive-beta.com

【この章の内容】


 この章で利用するパッケージを読み込みます。

# 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

 確率p1を、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))

f:id:anemptyarchive:20210206234423p:plain
ランダムウォーク

 試行回数$K$が大きくなるにつれて中心(平均値)から離れていくシミュレーション結果もあります。$2 \sqrt{K}$で補助線を引いています。どのくらいの確率で中心から離れるのかが気になった方は、二項分布を参照してください。ランダムウォークとの関係は書いてませんが、二項分布については「1.2.1:二項分布【『トピックモデル』の勉強ノート】 - からっぽのしょこ」でまとめています。

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))

f:id:anemptyarchive:20210206234453p:plain
連続値のランダムウォーク

 概ね$c \sqrt{K}$の範囲となるようです(どうやって求めるの?)。

参考文献

  • 花田政範・松浦壮『ゼロからできるMCMC』講談社,2020年.

おわりに

 既約性、詳細釣り合い条件についてはいい例がなかったので保留します。以前読んだ本に載ってたような、でもどの本かも覚えてない。

 さーんっ歩進んでっ2歩下がる~

【次の章の内容】

www.anarchive-beta.com