読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

R による棄却サンプリング法(おれが)入門

R 確率過程

はじめに;棄却サンプリング法とは

棄却サンプリング(rejection sampling)法は, 受容−棄却サンプリング(acceptance-rejection)法ともいう.

棄却サンプリング法を使うとなにができるか.

例えば, 密度関数 f(x) に従う乱数を発生させたいけど, 使っているソフトウェアにその乱数が用意されていないとき, 別の密度関数 g(x) に従う乱数から, f(x) に従う乱数を作ることができる.

乱数を発生させるためなら. 実際は逆関数法(逆関数法について知りたい方はググってください)のほうが使える場合が多いが, このエントリは最近人気の MCMC(Markov Chain Monte Calro; マルコフ連鎖モンテカルロ)を勉強するための練習, という感じ.

棄却サンプリング法のアルゴリズム

確率密度関数 g(x) に従う乱数を利用して, f(x) に従う乱数を作る.

g(x) に必要な仮定は, すべての点 x に対して

 \displaystyle f(x) \le c g(x)

が成立することである. ここで c は正の定数.

棄却サンプリング法のアルゴリズムは以下の通り.

  1. g(x) に従う乱数 x* を発生させる
  2. 区間 (0, 1) 上の一様乱数 u を発生させる
  3.  u \le f(x^\ast) / c g(x^\ast) ならば, x*f(x) からのランダム標本として受用(採用)し, そうでなければ棄却する(捨てる)
  4. これを好きなだけ(m 回)繰り返す

得られる x*f(x) に従う.

数値実験;カイ二乗分布からガンマ分布

めのこでやる

shape パラメータ = 1, scale パラメータ = 2 のガンマ分布を f(x), 自由度 1 のカイ二乗分布を g(x) としてガンマ分布に従う乱数を生成してみる.

 \displaystyle f(x) \le c g(x)

である必要がある.

まず c をめのこで決めてみる. c = 4 にしよう.

f:id:abrahamcow:20150112010932p:plain

 f(x) \le c g(x) が成立している.

f <- function(x){
  dgamma(x,shape=1,scale=2)
}

cg <- function(x){
 4 * dchisq(x,df=1)
}

curve(cg(x),0,4, lwd=2,lty=2, ylab="")
curve(f(x),add=TRUE,lwd=2)
legend("topright",c("f(x)", "cg(x)"),lty=1:2,lwd=2)

アルゴリズム通りに実行して, ヒストグラムで当てはまりを見る.

generator <- function(m){
  x <-rchisq(m,df=1)
  u = runif(m)
  flags <- u <= (f(x)/cg(x))
  acc <-sum(flags)
  return(x[flags])
}
m=10000
x <- generator(m)

hist(x, breaks="scott",freq=FALSE)
curve(f(x),lty=2,add=TRUE)
legend("topright","f(x)",lty=2,lwd=1)

f:id:abrahamcow:20150112021107p:plain

だいたい f(x)(ガンマ分布)に従う乱数が生成されたようだ.

ただし受容率(受容された割合)をみると, 約25%.

効率はあまりよくない.

> len <- length(x)
> len/m #受用率
[1] 0.2511

もう少しちゃんとやる

つぎに c をちゃんと計算してみる.

 f(x) \le c g(x) を知りたければ,  f(x)/g(x) \le c なる c を求めればいい.

ガンマ分布の密度関数は,

 \displaystyle x^{k-1} \frac{e^{-x/\theta}}{\Gamma(k)}\,\theta^k

なので,

shape パラメータ θ=1, scale パラメータ k =2 のガンマ分布は,

 \displaystyle f(x) =  x \exp(x),

カイ二乗分布の密度関数は,

 \displaystyle f(x;k)=\frac{(1/2)^{k/2}}{\Gamma(k/2)} x^{k/2 - 1} e^{-\frac{x}{2}}

 \Gamma(1/2)=\pi なので, 自由度 1 のカイ二乗分布

 g(x)=\frac{1}{\sqrt{2 \pi}} x^{-\frac{1}{2}}.

 \displaystyle f(x)/g(x) = \sqrt{2 \pi} x^\frac{3}{2}\exp(-x/2)

計算をかんたんにするために対数をとり, 微分して,

\displaystyle\frac{d}{dx} \log\{f(x)/g(x)\} = \log\{ \sqrt{2 \pi} x^\frac{3}{2}\exp(-x/2)\} \\
\displaystyle = \frac{3}{2} - \frac{1}{2}

これを = 0 とおいて x について解くと, x = 3.

 \displaystyle c = \sqrt{2 \pi} 3^\frac{3}{2}\exp(-3/2)

とすればより f(x) に近い cg(x) が得られるはずだ.

cg <- function(x){
  ( ( sqrt(2*pi)*3^(3/2) )* exp(-3/2) ) * dchisq(x,df=1)
}

curve(cg(x),0,4, lwd=2,lty=2, ylab="")
curve(f(x),add=TRUE,lwd=2)
legend("topright",c("f(x)", "cg(x)"),lty=1:2,lwd=2)

f:id:abrahamcow:20150112013915p:plain

ふたたび乱数を生成してみる.

x <- generator(m)
hist(x, breaks="scott",freq=FALSE)
curve(f(x),lty=2,add=TRUE)
legend("topright","f(x)",lty=2,lwd=1)

f:id:abrahamcow:20150112020919p:plain

> len <- length(x)
> len/m #受用率
[1] 0.3443

受容率は, 約34%. さきほどより上がった.
(それでも, やはりあまり効率がいいとはいえないと思うが)

組み込み乱数と較べる

R にはガンマ分布に従う乱数を生成する関数が用意されている.

##組み込みの乱数
set.seed(2014) #乱数の種を固定
x2 <- rgamma(len,shape=1,scale=2)

この x2 と棄却サンプリング法で生成したさきほどの x を較べる.

##組み込みの乱数
set.seed(2014)
x2 <- rgamma(len,shape=1,scale=2)

f:id:abrahamcow:20150112014702p:plain

x2 のほうが若干ガンマ分布に近いかなーという感じ.

Q-Q プロットを書いてみると,

#Q-Q plot
#nihongo()
par(mfrow=c(1,2))
Q <- qgamma(ppoints(len), shape=1, scale=2)
qqplot(Q, x, xlab="パーセンタイルの理論値")
abline(0,1,col="red2")
qqplot(Q, x2,xlab="パーセンタイルの理論値")
abline(0,1,col="red2")

f:id:abrahamcow:20150112014952p:plain

やはり x2 のほうが精度がいい.

棄却サンプリング法は g(x) の, f(x) を覆う近似の精度に依存するようだ.

参考文献

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

主に pp.159-161 を参照した.
※『計算統計学の方法』はていねいでいい本だと思うけど, ある程度数学好きな人向け. お買い求めの際はご注意を. 翻して言うと, 証明とか計算をちゃんとフォローしたい方にはおすすめです.

入門・演習 数理統計

入門・演習 数理統計

ガンマ分布, カイ二乗分布の密度関数とか, 覚えてなかったので参照した.

関連エントリ

R によるメトロポリス・ヘイスティング(MH)法入門 - 廿TT