メトロポリスヘイスティングによるサンプリング
木村奏
2025-01-01
はじめに
メトロポリスヘイスティングによるサンプリングの実装を行う。実装ではR6クラスを使ってオブジェクト指向的にコーディングをした。MCMCのサンプリング方法を学びたい方、およびオブジェクト指向をRで実装したい方が主な対象である。
\(y=3x+4\)に対して\(e \sim N(0,16)に従う\)ノイズを乗せたデータを生成し、傾き\(3\)が分からないという体のもと傾きをベイズ的に事後分布から推定する問題設定とした。尤度関数および事前分布はそれぞれ以下の通りである。
尤度関数 \[ p(D|a)=\Pi_{i=1}^{N}\frac{1}{\sqrt{2\pi\sigma^2}}\exp\Big(-\frac{(y_{i}-(ax_{i}+4))^2}{2\sigma^2}\Big) \] 事前分布 \[ p(a)=\frac{1}{\sqrt{2\pi×2000000}}\exp\Big(-\frac{a^2}{2×2000000}\Big) \]
パッケージの読み込み
R6クラスでのオブジェクト指向プログラミングを行うためここではパッケージ”R6”を読み込む。
library( R6 )データ生成
「はじめに」で述べた方法にもとづいて線形のデータを生成した。特にこだわりはないがデータは\(0-10\)の範囲で生成している。
############################ Data Generation ###################################
x <- seq( 0, 10, 0.1 )
y <- 3 * x + 4
y <- y + rnorm( length( x ), 0, 4 )
plot( x, y, pch = 19 )
事後分布の定義
尤度関数と事前分布の積を取る形で事後分布を求め、パラメータ\(a\)の関数としてdistというオブジェクトに格納した。
########################## Function Definition #################################
# likelihood function
like <- function( a ) {
sum <- 1
for( i in 1 : length( x ) ) {
sum <- sum * ( 1 / sqrt( 2 * pi * 4 ) * exp( - ( y[ i ] - a * x[ i ] - 4 )^2 / ( 2 * 4 ) ) )
}
return( sum )
}
# prior distribution
prio <- function( b ) {
1 / sqrt( 2 * pi * 2000000 ) * exp( - b ^ 2 / ( 2 * 2000000 ) )
}
# posterior distribution
dist <- function( c ) {
return( like( c ) * prio( c ) )
}MHクラスの定義
次にR6クラスを用いてメトロポリスヘイスティングによるサンプリングクラスを作成した。事後分布とメトロポリスヘイスティングの初期値をコンストラクタの入力とし、sample()というメソッドを実行するとサンプリングがsampleVecというフィールドに記録されていく仕組みになっている。またサンプル結果を簡単に可視化できると便利なので可視化用のvisualize()というメソッドを用意した。showSummary()とshowParInfo()はそれぞれ事後分布とMCMCサンプリングに関わる情報を表示するメソッドである。。
MH <-
R6Class( "MH",
public = list(
posterior = NA, # posterior function
N.samp = 100000, # number of sampling
N.int = 100, # interval
burn = 100, # burnout
init = NA, # initial value
sampleVec = NA, # sampling with N.int interval after burnout
sampleVecAll = NA, # All of sampling from scratch
plots = list(), # list to store plots
initialize = function( posterior, init ) {
self$posterior <- posterior
self$init <- init
},
# method to run MCMC sampling
sampling = function(){
self$sampleVecAll <- self$init
for ( j in 1 : self$N.samp ) {
if( j %% ( self$N.samp %/% 10 ) == 0 ) {
cat( paste( ( 10 * j %/% ( self$N.samp %/% 10 ) ), "% Completed\n", sep = "" ) )
}
# previous sampling
x1 <- self$sampleVecAll[ j ]
# next candidate sampling
x2 <- rnorm( 1, x1, 1 )
# store sampling to samp following sampling role
alpha <- self$posterior( x2 ) / self$posterior( x1 )
if ( alpha >= 1 ) {
self$sampleVecAll <- c( self$sampleVecAll, x2 )
} else {
p <- runif( 1, 0, 1 )
if ( p < alpha ) {
self$sampleVecAll <- c( self$sampleVecAll, x2 )
} else {
self$sampleVecAll <- c( self$sampleVecAll, x1 )
}
}
# store sampling to sampleVec if satisfies interval condition
if ( j >= self$burn ) {
if ( j %% self$N.int == 0 ) {
self$sampleVec <- c( self$sampleVec, self$sampleVecAll[ j + 1 ] )
}
}
}
self$sampleVec <- self$sampleVec[ - 1 ]
self$plots <- list(
function() plot( self$sampleVecAll, ylab = "Sampling", xlab = "Iteration", main = "all of sampling" ),
function() plot( self$sampleVec, ylab = "Sampling", xlab = "Iteration", main = "sampling after burnout with interval" ),
function() hist( self$sampleVec, main = "Histogram of Sampling", xlab = "Sampling" ),
function() acf( self$sampleVec, plot = TRUE, main = "sampling after burnout with interval")
)
},
# method to visualize sampling
visualize = function() {
if( length( self$sampleVecAll ) == 1 ){
stop( "please run sampling() beforehand" )
} else {
for( i in 1 : length( self$plots ) ) {
self$plots[[ i ]]()
if( i < length( self$plots ) ) {
readline( prompt = "次のプロットを出力するにはEnterを押してください" )
}
}
}
},
showSummary = function() {
cat( "summary of posterior distribution is as follows:\n" )
cat( "[1] mean: ", mean( self$sampleVec ), "\n" )
cat( "[2] variance: ", var( self$sampleVec ), "\n" )
},
showParInfo = function() {
cat( "basic information on MCMC sampling is as follows:\n" )
cat( paste( "[1] number of sampling: ", self$N.samp, "\n" ) )
cat( paste( "[2] interval of sampling: ", self$N.int, "\n" ) )
cat( paste( "[3] burnout: ", self$burn, "\n" ) )
}
)
)実行結果
以下では実際にこのMHクラスを用いてサンプリングをした結果を示す。
インスタントの作成とサンプリング
事後分布の関数とサンプリングの初期値を入れてインスタンスを作成する。sampling()メソッドを走らせることによりサンプリングが開始する。
MCMC <- MH$new( dist, 2)
MCMC$sampling()## 10% Completed
## 20% Completed
## 30% Completed
## 40% Completed
## 50% Completed
## 60% Completed
## 70% Completed
## 80% Completed
## 90% Completed
## 100% Completedパラメータ情報の表示
サンプリングをする際に設定したburnoutの数などのパラメータを表示するためのメソッドである。
MCMC$showParInfo()## basic information on MCMC sampling is as follows:
## [1] number of sampling: 1e+05
## [2] interval of sampling: 100
## [3] burnout: 100事後分布の平均と分散
次のサマリーにあるように大体\(a=3\)で推定されており正しく機能したことが分かる。
MCMC$showSummary## function() {
## cat( "summary of posterior distribution is as follows:\n" )
## cat( "[1] mean: ", mean( self$sampleVec ), "\n" )
## cat( "[2] variance: ", var( self$sampleVec ), "\n" )
## }
## <environment: 0x0000021b521d19a8>プロットによる確認
visualize()はサンプリングの結果を可視化するメソッドである。Enterを押すことで次のプロットに移行するようにしてある。一つ目がすべてのサンプリング、二つ目が100回目以降のサンプリングを100サンプルごとに表示した結果となっている。三つ目は二つ目をヒストグラムで表したもの、四つ目は自己相関係数である。
MCMC$visualize()
## 次のプロットを出力するにはEnterを押してください
## 次のプロットを出力するにはEnterを押してください
## 次のプロットを出力するにはEnterを押してください