library(dplyr)   # いろいろ
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(lattice) # 3次元プロット
library(coda)    # Rhatの計算

等高線図(contour plot)を描いてみる

まずは等高線図のデータを格納するdata.frameを作成し、contour.dfと名付ける。このdata.frameにはポテンシャルエネルギーの母数(\(\theta\))thetaと運動エネルギーの母数(\(p\))pが入る。

contour.x <- seq(-10, 10, by = 0.2) 
contour.y <- seq(-10, 10, by = 0.2) 
contour.df <- expand.grid(contour.x, contour.y) # グリッド生成
colnames(contour.df) <- c("theta", "p")         # 列名の変更

グリッドが作成されたら、格納されているthetaとpを用いてハミルトニアンを計算する。ハミルトニアンは

\[H(\theta, p) = U(\theta) + K(p)\]

である。

ポテンシャルエネルギー

ポテンシャルエネルギ\(U(\theta)\)は\(U(\theta) = mgh(\theta)\)であるが、質量(\(m\))と重力加速度(\(g\))は定数1にすると無視できるので、\(h(\theta)\)のみを用いる。この\(h(\theta)\)は、事後分布を対数変換したうえで-1を乗じたものである。ここでは\(x = 2, \sigma = 1\)の正規分布を事後分布とし、母数\(\theta\)(普通なら\(\mu\)だが…)を推定したいとする。

\[p(\theta | 2) = \frac{1}{\sqrt{2 \pi}\sigma}\exp(-\frac{(2 - \theta)^2}{2 \sigma^2})\]

\(\sigma = 1\)だし、定数も全部取りぬくとカーネルは

\[\exp(-\frac{(2 - \theta)^2}{2})\]

ここに対数と-1をかけると

\[log(\exp(-\frac{(2 - \theta)^2}{2})) \times -1 = \frac{(2 - \theta)^2}{2}\]

になる。

運動エネルギー

運動エネルギー\(K(p)\)は\(K(p) = \frac{1}{2m}p^2\)であるが、ここでも質量もまた定数1とし、\(0.5 p^2\)のみを使う。

上のポテンシャルエネルギーと運動エネルギーの和がハミルトニアンである。このハミルトニアンは等高線図において「高さ」を意味することになる。

式も分かったし、早速、ハミルトニアンを計算するcontour.h関数を作成する。

# 高さを計算するcontour.h関数を作成
contour.h <- function(theta, p){
    PE <- 0.5 * (2 - theta)^2  # Potential Energy
    KE <- 0.5 * p^2            # Kinetic Energy
    return(PE + KE)            # Hamiltonian: P.E. * K.E.
}

先ほど作成したcontour.dfに新しい列hを追加し、contour.h関数を使ってハミルトニアンを計算し、格納する。

# contour.dfにh列を生成し、contour.h関数を用いて
# 等高線図における高さ(=ハミルトニアン)を格納する
contour.df <- contour.df %>% mutate(h = contour.h(theta, p))

等高線図と3次元プロットを確認する。

# 等高線図の表示
contour(x = seq(-10, 10, 0.2), y = seq(-10, 10, 0.2), 
        z = matrix(contour.df$h, nrow = 101, byrow = TRUE),
        xlab = "p", ylab = "theta")

# 3Dプロット。Latticeの使い方が分からないので
# ネットで適当にコードを持ってきました。
wireframe(h ~ theta * p, data = contour.df,
  xlab = "theta", ylab = "p", zlab = "Hamiltonian",
  main = "Surface Plot",
  drape = TRUE,
  colorkey = TRUE,
  screen = list(z = -60, x = -60)
)

もうちょっと気持ち悪い形の等高線図が欲しかったが、仕方ない…

本格的にHMCをやってみる

# ハミルトニアンを計算する関数
lf_hamil <- function(theta, p){
    PE <- 0.5 * (2 - theta)^2  # Potential Energy
    KE <- 0.5 * p^2            # Kinetic Energy
    return(PE + KE)            # Hamiltonian: P.E. * K.E.
}

# ハミルトニアン式を微分
lf_devi_theta  <- function(theta){
    return(theta - 2)
}

# pを0.5ステップ移動
lf_halfstep_p  <- function(theta, p, e = 0.01){
    return(p - 0.5 * e * lf_devi_theta(theta))
}

# thetaを1ステップ移動
lf_1step_theta <- function(theta, p, e = 0.01){
    return(theta + e * p)
}

# L回移動
lf_1step <- function(theta, p, e = 0.01, L = 100, gap = 1){
    return.df <- data.frame(matrix(rep(NA, (L+1) * 4), 
                            ncol = 4))
    colnames(return.df) <- c("id", "theta", "p", "h")
    return.df[1, ] <- c(1, theta, p, lf_hamil(theta, p))
        
    for(i in 2:(L+1)){
        p     <- lf_halfstep_p(theta, p, e)
        theta <- lf_1step_theta(theta, p, e)
        p     <- lf_halfstep_p(theta, p, e)
        return.df[i, ] <- c(i, theta, p, lf_hamil(theta, p))   
    }

    class(return.df) <- append("Lstep", class(return.df))

    return(return.df)
}

# 1ステップ×L回の結果を図示
plot.Lstep <- function(df){
    contour(seq(-10, 10, 0.2), seq(-10, 10, 0.2), 
            matrix(contour.df$h, nrow = 101, byrow = TRUE),
         xlab = "p", ylab = "theta")
    lines(x = df$p, y = df$theta)
    points(x = df$p, y = df$theta, pch = 19,
           cex = 0.5)
    text(x = df$p[1], y = df$theta[1], 
         labels = "Start", pos = 1)
    text(x = df$p[nrow(df)], y = df$theta[nrow(df)],
         labels = "End", pos = 1)
}

まず、テストとして\(theta = 4, p = -5, \varepsilon = 0.05, L = 100\)の設定で移動させ、その結果を見る。また、その移動の軌跡を図で確認する。

test_step <- lf_1step(theta = 4, p = -5, e = 0.05, L = 100)

head(test_step, 20)
##    id       theta         p        h
## 1   1  4.00000000 -5.000000 14.50000
## 2   2  3.74750000 -5.093687 14.49970
## 3   3  3.49063125 -5.174641 14.49944
## 4   4  3.23003592 -5.242657 14.49922
## 5   5  2.96636550 -5.297567 14.49904
## 6   6  2.70027917 -5.339234 14.49890
## 7   7  2.43244214 -5.367552 14.49881
## 8   8  2.16352401 -5.382451 14.49876
## 9   9  1.89419706 -5.383894 14.49875
## 10 10  1.62513463 -5.371877 14.49879
## 11 11  1.35700935 -5.346431 14.49888
## 12 12  1.09049155 -5.307618 14.49901
## 13 13  0.82624753 -5.255537 14.49918
## 14 14  0.56493788 -5.190316 14.49939
## 15 15  0.30721589 -5.112120 14.49965
## 16 16  0.05372586 -5.021144 14.49993
## 17 17 -0.19489848 -4.917614 14.50026
## 18 18 -0.43803558 -4.801791 14.50061
## 19 19 -0.67507759 -4.673963 14.50099
## 20 20 -0.90543191 -4.534451 14.50139
plot(test_step)

head(test_step, 20)の結果をみると\(\theta\)と\(p\)が動いてもハミルトニアン(\(h\))はほぼ一定だということが分かる。これは\(\varepsilon\)の大きさの関係がある。

前回の勉強会ではeとLを変えたらどうなるかで議論があったが、ここで実際にやってみる。やってみる設定は

である。pとthetaの初期値は比較のためにそのまま使う。

順次にやってみよう。

theta <-  4 # thetaの初期値
p     <- -5 # pの初期値

lf_test1 <- lf_1step(theta, p, e = 0.01, L = 100) # 設定1
lf_test2 <- lf_1step(theta, p, e = 0.1,  L = 100) # 設定2
lf_test3 <- lf_1step(theta, p, e = 0.5,  L =  10) # 設定3
lf_test4 <- lf_1step(theta, p, e = 0.75, L =   5) # 設定4

head(lf_test1, 10) # 設定1の結果
##    id    theta         p        h
## 1   1 4.000000 -5.000000 14.50000
## 2   2 3.949900 -5.019749 14.50000
## 3   3 3.899605 -5.038997 14.50000
## 4   4 3.849120 -5.057741 14.49999
## 5   5 3.798450 -5.075979 14.49999
## 6   6 3.747600 -5.093709 14.49999
## 7   7 3.696576 -5.110930 14.49999
## 8   8 3.645382 -5.127639 14.49998
## 9   9 3.594023 -5.143836 14.49998
## 10 10 3.542505 -5.159519 14.49998
plot(lf_test1)     # 設定1の図示

head(lf_test2, 10) # 設定2の結果
##    id      theta         p        h
## 1   1  4.0000000 -5.000000 14.50000
## 2   2  3.4900000 -5.174500 14.49778
## 3   3  2.9651000 -5.297255 14.49616
## 4   4  2.4305490 -5.367037 14.49523
## 5   5  1.8916925 -5.383150 14.49501
## 6   6  1.3539191 -5.345430 14.49552
## 7   7  0.8226065 -5.254256 14.49673
## 8   8  0.3030678 -5.110540 14.49860
## 9   9 -0.1995015 -4.915718 14.50105
## 10 10 -0.6800759 -4.671740 14.50398
plot(lf_test2)     # 設定2の図示

head(lf_test3, 10) # 設定3の結果
##    id     theta          p        h
## 1   1  4.000000 -5.0000000 14.50000
## 2   2  1.250000 -5.3125000 14.39258
## 3   3 -1.312500 -4.2968750 14.71790
## 4   4 -3.046875 -2.2070312 15.17097
## 5   5 -3.519531  0.4345703 15.32704
## 6   6 -2.612305  2.9675293 15.03979
## 7   7 -0.552002  4.7586060 14.57852
## 8   8  2.146301  5.3600311 14.37567
## 9   9  4.808029  4.6214485 14.62141
## 10 10  6.767750  2.7275038 15.08536
plot(lf_test3)     # 設定3の図示

head(lf_test4, 10) # 設定4の結果
##   id      theta         p        h
## 1  1  4.0000000 -5.000000 14.50000
## 2  2 -0.3125000 -4.882812 14.59476
## 3  3 -3.3242188 -2.019043 16.21192
## 4  4 -3.3410645  1.980438 16.22455
## 5  5 -0.3535614  4.865923 14.60823
## 6  6  3.9578199  5.014326 14.48826
plot(lf_test4)     # 設定4の図示

つまり、等高線にどれだけ滑らかに沿って行くかを決めるのがeで、何回行くかを決めるのがLである。eが非常に小さければ、同じハミルトニアンを維持することになるが、大きくなるとズレも大きくなる。

\(e\)の大きさとハミルトニアンの変動も見てみる。設定1は\(e = 0.01\)、設定2は\(e = 0.1\)であり、そのハミルトニアン(\(h\))の分散は

var(lf_test1$h)
## [1] 1.140514e-09
var(lf_test2$h)
## [1] 0.0001729569

設定1(\(e = 0.01\))の\(h\)の分散は約\(1.14^{-9}\)、設定2(\(e = 0.1\))は約\(1.73^{-4}\)であり、\(e\)が小さい方がよりハミルトニアンが維持されることが分かる。

今までのステップは複数回繰り返す

run_hmc <- function(theta = NULL, e = 0.05, 
                    L = 100, iter = 100){
    # HMCの実行
    # 引　数: theta - thetaの初期値
    #     　  e     - 滑らかさ(?)
    #   　    L     - Leapの回数
    #    　   iter  - 繰り返し回数
    # 参　考: thetaがNULLなら、ランダムに決定
    # 返り値: List型
    
    # 関数の引数を保存
    init.set <- c(theta, e, L, iter)
    
    # thetaが指定されていないなら初期値はUnif(-10, 10)から
    if(is.null(theta)) theta <- runif(1, min = -10, max = 10)
    # pの初期値
    p <- rnorm(1, mean = 0, sd = 1)
    
    # 1ステップ×L回後のthetaとp値を格納するdata.frame
    result.df <- data.frame(matrix(rep(NA, iter * 8), 
                                   ncol = 8))
    colnames(result.df) <- c("id", "theta", "p", 
                             "prev_theta",
                             "prev_H", "curr_H", "r",
                             "Accept")
    
    # result.dfの1行目に初期値を格納
    result.df[1, ] <- c(1, theta, p,        # ID, theta, p
                        NA,                 # previous theta
                        NA,                 # previous Hamiltonian
                        lf_hamil(theta, p), # Current Hamiltonian
                        NA,                 # r
                        1)                  # Accepted?
    
    # 結果を格納するリストを作成
    result.list <- list()
    
    i       <- 1 # Iteration index
    i.count <- 1 # Iteration index(Accepted only)
    
    while(i.count <= iter){

        p      <- rnorm(1, mean = 0, sd = 1) # pを生成
        prev_p <- p                          # Leapする前のpを保存
        
        # t-1でthetaが受容されたら前回のLeap後のthetaが今回のthetaになる
        # コード上prev_thetaとthetaは同じ値なので、区分する必要はないが、
        # 理解のためにわざと別けておく
        if(result.df$Accept[i] == 1){
            prev_theta <- result.df$theta[i]
            theta      <- result.df$theta[i]
        }else{ # 受容されなかったら、受容された最後のthetaが今回のthetaになる
            prev_theta <- result.df$prev_theta[i]
            theta      <- result.df$prev_theta[i]
        }
        
        # prev_thetaとprev_pのハミルトニアン
        prev_H     <- lf_hamil(prev_theta, prev_p)
        
        # ジャンプ! ジャンプ!
        onestep <- lf_1step(theta, p, e, L)
        
        # L回ステップの軌跡を保存しておく (いつか使えるかも)
        result.list[[i]] <- onestep
        
        # L回ステップ後のthetaとpを用いてハミルトニアンを計算する
        curr_H <- lf_hamil(onestep$theta[L+1],
                           onestep$p[L+1])
        
        # 受容確率を計算
        r      <- exp(prev_H - curr_H)
        
        # rが1より小さく、ランダム値がrより大きかったら受容されず、
        # それ以外は受容される
        if(r < 1 & runif(1, 0, 1) > r){
            accept  <- 0           # Not accepted...
            i.count <- i.count     # No increasing the counter
        }else{
            accept  <- 1           # Accepted!
            i.count <- i.count + 1 # Increasing the counter
        }
        
        # result.list[[1]]に結果を保存
        # result.list[[2]]移行は全ての軌跡が保存されるが、
        # ここではL回ステップ後の結果のみが格納される。
        result.df[i+1, ] <- c(i+1,                 # ID
                              onestep$theta[L+1],  # ステップ後のtheta
                              onestep$p[L+1],      # ステップ後のp
                              prev_theta,          # 前回のtheta
                              prev_H, curr_H,      # 前回と今回のハミルトニアン
                              r, accept)           # 受容確率と受容結果
        
        i <- i + 1
    }
    
    result.list[["overall"]] <- result.df
    result.list[["setting"]] <- init.set
    
    # 単にクラスとメソッドが使ってみたい
    class(result.list) <- append(class(result.list), "HMC")
    
    return(result.list)
}
plot.HMC <- function(result_list, alpha = 1,
                     only_accept = TRUE,
                     hist_break = 30,
                     hist_adj = 1){
    # hmc_plot_simple
    # HMC結果(Leapは省略)の図示
    # 引数: result_list - 結果オブジェクト(List型)
    #     : only_accept - 受容された結果のみを表示
    #     : alpha       - Contour Plotの線の透明度
    #     : hist_breaks - Histgoramのbinの数
    #     : hist_adj    - Density Lineの調整
    
    # MAP, EAP, MED, 受容率を計算
    result_summary <- summary(result_list)
    
    # ACF & ESS
    acf_ess        <- ESS(result_list)
    
    # HMCのセッテイングを読み込む
    init.set       <- result_list[["setting"]]
    
    # 受容された結果のみ示すならAcceptが1のケースのみを抽出
    if(only_accept == TRUE){
        df    <- result_list[["overall"]] %>% filter(Accept == 1)
        df$id <- 1:nrow(df)
    }else{
        df    <- result_list[["overall"]]
    }
    
    
    par(mfrow = c(2, 2), oma = c(0, 0, 2, 0)) # 等高線図、thetaのヒストグラム、trace plot
    # 1. 等高線図
    contour(seq(-10, 10, 0.2), seq(-10, 10, 0.2), 
            matrix(contour.df$h, nrow = 101, byrow = TRUE),
            xlab = "p", ylab = "theta", 
            main = "Contour Plot")
    lines(x = df$p, y = df$theta,
          col = rgb(0, 0, 0, alpha))
    points(x = df$p, y = df$theta,
           pch = 19, cex = 0.5, col = rgb(0, 0, 0, alpha))
    
    # 2. ヒストグラム
    hist(df$theta, breaks = hist_break,
         probability = TRUE,
         xlab = "theta", ylab = "Density",
         main = "Histogram of theta")
    lines(density(df$theta, adjust = hist_adj), 
          col="red", lwd = 2)
    legend("topright", legend = c(paste("MAP:", result_summary$MAP, " "), 
                                  paste("EAP:", result_summary$EAP, " "), 
                                  paste("MED:", result_summary$MED, " ")),
           bty = "n", cex = 1, pt.cex = 1)
    
    # 3. Trace plot
    plot(x = df$id, y = df$theta,
         type = "l", xlab = "ID", ylab = "theta", 
         main = "Trace Plot")
    legend("topright", bty = "n", cex = 1, pt.cex = 1,
           legend = c(paste("Accept Ratio:", 
                            round(result_summary$Accept_Ratio, 5))))
    
    # 4. Autocorrelation
    plot(x = 0:(length(acf_ess$ACF)-1), y = abs(acf_ess$ACF), 
         type = "h", main = "Autocorrelation", lwd = 3,
         col = "skyblue", xlab = "Lag", ylab = "ACF")
    abline(h = 0)
    legend("topright", bty = "n", cex = 1, pt.cex = 1,
           legend = c(paste("ESS:", round(acf_ess$ESS, 5))))
    
    # Title
    mtext(paste("theta = ", init.set[1], "; e =", init.set[2],
                "; L = ", init.set[3], "; iter =", init.set[4]), outer=TRUE,
          cex = 1.5)
    par(mfrow = c(1, 1), oma = c(0, 0, 0, 0))
}

HMC結果の要約を返す関数summary.HMCを作成する。run_HMC()から返されるオブジェクトはHMCクラスなので、実際にはsummary(オブジェクト名)で作動する。

返される結果はEAP, MAP, MED, 受容確率であり、List型である。

summary.HMC <- function(result_list, burn = NULL){
    # summary.HMC
    # HMC結果の統計量
    # 引　数: result_list - 結果オブジェクト(List型)
    #   　    burn        - burn-in数
    # 返り値: List型
    # 参　考: burnがNAなら最初の10%を捨てる
    
    # 受容されたケースのみを抽出
    df <- filter(result_list[["overall"]], Accept == 1)
    
    # Burn-in期間が設定されてないなら、最初の10%を切る
    if(is.null(burn)){
        burn <- round(nrow(df) * 0.1, 0)
    }
    
    result_theta <- df$theta[-(1:burn)]
        
    # EAP, MED, MAPを計算
    EAP <- round(mean(result_theta), 5)
    MED <- round(median(result_theta), 5)
    MAP <- round(density(result_theta)$x[which.max(density(result_theta)$y)], 5)
    
    # 受容確率
    AR  <- mean(result_list[["overall"]]$Accept)
    
    # 結果を返す
    return(list(EAP = EAP,             # EAP(平均値)
                MED = MED,             # MED(中央値)
                MAP = MAP,             # MAP(最頻値)
                Accept_Ratio = AR))  # 受容確率
}

続いて、ESSを計算するESS.HMC関数を作成する。これもまたESS(オブジェクト名)で走らせる。

ESSの計算法はKruschkeのp.184を参照した。

ESS <- function(x) { 
    UseMethod("ESS", x) 
}

ESS.HMC <- function(result_list, max.lag = 20){
    # ESS.HMC
    # ACFとESSを計算
    # 引数: result_list - 結果オブジェクト(List型)
    #       max.lag     - 最大ラグ
    
    # 受容されたケースのみを抽出
    theta     <- filter(result_list[["overall"]], Accept == 1)
    theta     <- theta$theta
    theta_len <- length(theta)
    
    acf.vector <- as.vector(acf(theta, plot = FALSE, 
                                lag.max = theta_len)$acf)
    result.vec <- abs(acf.vector)
    
    # ACF < 0.05 になるlagポイントを特定する
    ESS.vec <- c()
    for(i in 1:theta_len){
        if(result.vec[i] < 0.05){
            break
        }
        ESS.vec[i] <- result.vec[i]
    }
    
    # ESS
    ESS <- theta_len/(1 + 2 * sum(ESS.vec))
    
    # 結果を返す
    return(list(ACF = acf.vector[1:(max.lag+1)], # ACF (Vector)
                ESS = ESS))                      # ESS (Scala)
}

早速HMCをやってみる。まずは母数(\(\theta\))の初期値は10、 eは10、Lは100とし、これを100回繰り返す。

test_hmc <- run_hmc(theta = 10, e = 0.01, 
                    L = 100, iter = 100)

要約統計量(点推定値と受容確率)とプロットを表示する。

(要約統計量はプロット内でも表示される。)

summary(test_hmc)
## $EAP
## [1] 1.80174
## 
## $MED
## [1] 1.91037
## 
## $MAP
## [1] 2.30224
## 
## $Accept_Ratio
## [1] 1
plot(test_hmc, alpha = 0.25)

問題なく動くようなので、いろいろとパラメータをいじめてみる。

eとLが与える影響

e = 0.1 e = 0.5 e = 0.9
L = 5 設定1 設定2 設定3
L = 50 設定4 設定5 設定6
L = 100 設定7 設定8 設定9 
eL_test1 <- run_hmc(theta = 10, e = 0.1, # 設定1
                    L = 5, iter = 3000)
eL_test2 <- run_hmc(theta = 10, e = 0.5, # 設定2
                    L = 5, iter = 3000)
eL_test3 <- run_hmc(theta = 10, e = 0.9, # 設定3
                    L = 5, iter = 3000)
eL_test4 <- run_hmc(theta = 10, e = 0.1, # 設定4
                    L = 50, iter = 3000)
eL_test5 <- run_hmc(theta = 10, e = 0.5, # 設定5
                    L = 50, iter = 3000)
eL_test6 <- run_hmc(theta = 10, e = 0.9, # 設定6
                    L = 50, iter = 3000)
eL_test7 <- run_hmc(theta = 10, e = 0.1, # 設定7
                    L = 100, iter = 3000)
eL_test8 <- run_hmc(theta = 10, e = 0.5, # 設定8
                    L = 100, iter = 3000)
eL_test9 <- run_hmc(theta = 10, e = 0.9, # 設定9
                    L = 100, iter = 3000)
plot(eL_test1, alpha = 0.25)

plot(eL_test2, alpha = 0.25)

plot(eL_test3, alpha = 0.25)

plot(eL_test4, alpha = 0.25)

plot(eL_test5, alpha = 0.25)

plot(eL_test6, alpha = 0.25)

plot(eL_test7, alpha = 0.25)

plot(eL_test8, alpha = 0.25)

plot(eL_test9, alpha = 0.25)

ESSの比較

e = 0.1 e = 0.5 e = 0.9
L = 5 71.09522 53.28369 333.66667
L = 50 269.13347 6.42928 67.98044
L = 100 78.81037 23.53219 189.72671

ESSが大きい方が良いとは言われているが、どれくらいがいいかはモデルによる。

eが大きく、Lが小さくなるとESSが高くなる傾向があるが、これは受容確率に悪影響を与える。詳細はこれから

受容確率の比較

e = 0.1 e = 0.5 e = 0.9
L = 5 1.00000 0.98901 0.92408
L = 50 1.00000 0.99600 0.97003
L = 100 0.99700 0.99301 0.93207

結果を見ると受容確率にはeとL両方が影響を与えると考えられる。その中でeの影響力が高く、Lはそこまで大きくはない。全体的な傾向としては

  • eが小さいほど受容確率は上がる。
  • Lが小さいほど受容確率は上がる。(ただし、二乗の関係があるかも)

受容確率は1に近いほど良いと言われているし、ならばeとLは小さい方がいいのか。続いて、収束の速度も一緒に考えて見よう。

収束診断

収束診断のために最も一般的に用いられる指標として\(\hat{R}\)がある。これはチェーン間の分散とチェーン内の分散の比である。これまで使ったコードはチェーン数が1なので、\(\hat{R}\)を求めることは出来ない。したがって、先ほど作ったrun_hmc関数を使って複数のチェーンでHMCを行う関数、run_multiple_hmcを作成する。この関数は\(\hat{R}\)を見るためのものなので、thetaだけを保存し、全ての情報を捨てる。

ここではcodaパッケージを使う。MCMCから得られた値をmcmcあるいはmcmc.listクラスで保存しておくと色々使い道がある。

multiChain <- function(x, chains){
    UseMethod("multiChain", x)
}

multiChain.HMC <- function(HMC.list, chains = 2){
    
    # multiChain.HMC
    # HMCを複数回繰り返す
    # 依　存: dplyr, coda
    # 引　数: HMC.list - HMCオブジェクト
    # 　      chain    - chain数
    # 返り値: mcmc.list
    
    # HMCオブジェクトから初期値設定を読み込む
    init.set <- HMC.list[["setting"]]
    
    # 結果を格納するリストを作成
    result_list <- list()
    
    # 結果リスト1は既にあるからそのまま使う
    temp.df    <- HMC.list$overall %>% filter(Accept == 1)
    
    result_list[[1]] <- mcmc(temp.df$theta) # codaパッケージのmcmcオブジェクト
    
    # (指定したチェーン数 - 1)回だけ反復する。
    for(i in 2:chains){
        temp.list  <- run_hmc(init.set[[1]], init.set[[2]], 
                              init.set[[3]], init.set[[4]])
        temp.df    <- temp.list$overall
        temp.df    <- temp.df %>% filter(Accept == 1)
        
        result_list[[i]] <- mcmc(temp.df$theta) # codaパッケージのmcmcオブジェクト
    }
    
    return(mcmc.list(result_list)) # codaパッケージのmcmc.listオブジェクト
}

上記の９つの設定で10-chainのHMCを実行する。

eL_testオブジェクト内にはHMCの初期値も格納されているため、それを読み込んで自動的に同じ設定で指定したチェーン数-1回反復する。

# 同じ設定で9回やる(1回目は既にやったので)
mcTest1 <- multiChain(eL_test1, chains = 10)
mcTest2 <- multiChain(eL_test2, chains = 10)
mcTest3 <- multiChain(eL_test3, chains = 10)
mcTest4 <- multiChain(eL_test4, chains = 10)
mcTest5 <- multiChain(eL_test5, chains = 10)
mcTest6 <- multiChain(eL_test6, chains = 10)
mcTest7 <- multiChain(eL_test7, chains = 10)
mcTest8 <- multiChain(eL_test8, chains = 10)
mcTest9 <- multiChain(eL_test9, chains = 10)

codaパッケージ内のgelman.plotを使えば\(\hat{R}\)のプロットが表示される。

この(\(\hat{R}\))はまた使うので、保存しておこう。

# Rhat値を格納する空Listを作成
Rhat_list <- list()

# プロットの表示と同時にRhatの値は保存しておく
Rhat_list[[1]] <- gelman.plot(mcTest1, bin.width = 1, max.bins = 3000)$shrink[,,1]

Rhat_list[[2]] <- gelman.plot(mcTest2, bin.width = 1, max.bins = 3000)$shrink[,,1]

Rhat_list[[3]] <- gelman.plot(mcTest3, bin.width = 1, max.bins = 3000)$shrink[,,1]

Rhat_list[[4]] <- gelman.plot(mcTest4, bin.width = 1, max.bins = 3000)$shrink[,,1]

Rhat_list[[5]] <- gelman.plot(mcTest5, bin.width = 1, max.bins = 3000)$shrink[,,1]

Rhat_list[[6]] <- gelman.plot(mcTest6, bin.width = 1, max.bins = 3000)$shrink[,,1]

Rhat_list[[7]] <- gelman.plot(mcTest7, bin.width = 1, max.bins = 3000)$shrink[,,1]

Rhat_list[[8]] <- gelman.plot(mcTest8, bin.width = 1, max.bins = 3000)$shrink[,,1]

Rhat_list[[9]] <- gelman.plot(mcTest9, bin.width = 1, max.bins = 3000)$shrink[,,1]

一つずつ見ていくと比較しづらい。

比較のために同じプロットで表示させてみよう。

# 9つの設定のHMCから得られたRhatを同時にプロット
plot(x = 51:3002, y = Rhat_list[[1]], 
     type = "l", yaxt = "n", ylim = c(1, 2), col = 1,
     xlab = "Iteration", ylab = "Rhat")
for(i in 2:length(Rhat_list)){
    lines(x = 51:3002, y = Rhat_list[[i]], col = i)
}
axis(1, at=seq(1, 2, by = 0.1), labels = FALSE)
text(y = seq(1, 2, by = 0.1), 0, labels = seq(1, 2, by = 0.1), srt = 0, pos = 2, xpd = TRUE)
abline(h = 1.1, lty = 2)
legend("topright", lty = rep(1, 9), lwd = rep(1, 9), col = 1:9, 
       legend = c("e 0.1; L 5",   "e 0.5; L 5",   "e 0.9; L 5",
                  "e 0.1; L 50",  "e 0.5; L 50",  "e 0.9; L 50",
                  "e 0.1; L 100", "e 0.5; L 100", "e 0.9; L 100"))

ほとんどの設定で収束の速度はかなり早いが、

  • 設定1: e = 0.1; L = 5
  • 設定5: e = 0.5; L = 50
  • 設定8: e = 0.5; L = 100

の設定では比較的に収束が遅い。とりわけ、設定5は2000回まで収束していない。

この結果からは体系的な関係は見られない。多分、eとLが独立的に収束速度に影響を与えるというよりは、組み合わせ(交互作用)で影響を与えているのではないかと考えられる。

eとLの設定で収束の時間が変わるならM-Hと比べ、何がいいのか。