「Rによるモンテカルロ法入門」読書ノート：2章逆変換・一般変換法

読書ノートアーカイブ：http://yagays.github.com/blog/2012/10/20/archive-introducing-monte-carlo-methods-with-r/

内容：一様分布を元にして様々な確率分布に従う乱数を生成する

統計解析ソフトであるRには，様々な分布に対応した組み込み関数が用意されている．本章では，そういった組み込み関数を使わずに，一様分布から生成される乱数を逆関数で変換することで，他の確率分布の乱数を表現する．

本書で扱う乱数とは，完全なランダム性を持つ乱数ではなく擬似乱数である．擬似乱数はset.seed()関数で設定した値を種として乱数を生成するため，どのような環境においてもset.seed()関数で同じ値を使うことで乱数を再現することができる．

今回は，一様乱数を元にして指数乱数を生成し，指数乱数からガンマ分布やベータ分布の乱数へと変換していく．これらの確率分布は全てRの関数で用意されているものなので，一様乱数から生成した乱数とRの関数から生成した乱数を比較することによって，乱数の生成がうまくいっているかどうかを判断する．以下のコードで描写したヒストグラムはすべて，左側が一様分布から生成した乱数，右側がRの関数を用いて生成した乱数(N=10⁴ )．赤い曲線はどちらもRの関数を用いて分布を示したもの．

例 2.1 一様乱数から指数乱数を作成する

これはテキストにある通りのコードと作図．mcsmパッケージのdemo(Chapter.2)にも同様のコードがある．

以下のコードでは $-\log{(1-u)}$ をそのままコードに落とし込んでいる．テキストにもある通り， $U \sim \mathcal{U}_{[0,1]}$ ならば0から1の間で一様に等しい確率分布なのだがら， $U$ も $1-U$ も一様分布になる．

par(mfrow=c(1,2))
Nsim <- 10^4
U <- runif(Nsim)
X <- -log(1-U)
Y <- rexp(Nsim)
par(mfrow=c(1,2))
hist(X, freq=F, main="Exp from Uniform", nclass=50)
curve(dexp(x), add=T, col="red", lwd=2)
hist(Y, freq=F, main="Exp from R", nclass=50)
curve(dexp(x), add=T, col="red", lwd=2)

練習問題 2.2 逆変換法を用いてロジスティック分布とコーシー分布の乱数を生成する

逆変換法を用いてロジスティック分布とコーシー分布の乱数を生成する．変換で求まる関数はいわゆる逆累積分布関数というもの．

a.ロジスティック分布

$u = \frac{1}{1+e^{-\frac{x-\mu}{\beta}}}$ を変形して， $x = -\beta \log{\frac{1-u}{u}} + \mu$

以下のコードは $\mu=0$ ， $\beta=1$ の場合．

Nsim <- 10^4
U <- runif(Nsim)
par(mfrow=c(1,2))
X <- -log((1-U)/U)
Y <- rlogis(Nsim, 0, 1)
hist(X,freq=F,main="Logis from Uniform")
curve(dlogis(x), add=T, col="red", lwd=2)
hist(Y, freq=F, main="Logis from R")
curve(dlogis(x), add=T, col="red", lwd=2)

b.コーシー分布

$u = \frac{1}{2} + \frac{1}{\pi} \arctan{\frac{x-\mu}{\sigma}}$ を変形して， $x = \mu + \sigma \tan{\pi(u-\frac{1}{2})}$

以下のコードは $\mu=0$ ， $\sigma=1$ の場合．

Nsim <- 10^4
U <- runif(Nsim)
x <- tan(pi*(U-0.5))
Y <- rcauchy(Nsim)
hist(X, freq=F, main="Cauchy from Uniform", xlim=c(-10,10))
curve(dcauchy(x), add=T, col="red", lwd=2)
hist(Y[Y>=-10 & Y<=10], freq=F, main="Cauchy from R", xlim=c(-10,10))
curve(dcauchy(x), add=T, col="red", lwd=2)
# cauchyは値が両端に飛んでヒストグラムが綺麗に書けないので[-10,10]で整えている

練習問題 2.12 指数分布からガンマ分布とベータ分布の乱数を生成する

a.ガンマ乱数

以下の作図では， $\beta$ の値を固定して $a$ の値を1,2,5,9と変化させたときの分布の変化を見ている．

Nsim <- 10^4
beta <- 2.0
par(mfrow=c(4,2))
for(a in c(1,2,5,9)){
  U <- runif(a*Nsim)
  m <- matrix(U, nrow=a)
  X <- beta * apply(-log(m), 2, sum)
  Y <- rgamma(Nsim, shape=a, scale=beta)
  hist(X, freq=F, main=paste("Gamma from Uniform:","a=",a," beta=",beta), nclass=50)
  curve(dgamma(x, shape=a, scale=beta), add=T, col="red", lwd=2)
  hist(Y, freq=F, main=paste("Gamma from Uniform:","a=",a," beta=",beta), nclass=50)
  curve(dgamma(x, shape=a, scale=beta), add=T, col="red", lwd=2)
}

a.ベータ乱数

この手法では $(a \in \mathbb{N^*})$ ， $\mathbb{N^*} = 1,2, \dots$ という制約があるため， $a=0.1$ ， $b=0.1$ のようなベータ分布は作ることが出来ない．以下はa=1,b=1，a=2,b=3，a=8,b=4の場合．

# a=1,b=1の場合
par(mfrow=c(3,2))
Nsim <- 10^4
a <- 1
b <- 1
U <- runif((a+b)*Nsim)
m <- matrix(U,nrow=(a+b))
X <- -log(m[a,]) / apply(-log(m),2,sum) # a>=2の場合はapply(-log(m[1:a,]),2,sum) / apply(-log(m),2,sum)
Y <- rbeta(Nsim,a,b)
hist(X, freq=F, main=paste("Beta from Uniform:","a=",a,"b=",b), nclass=50, ylim=c(0,3))
curve(dbeta(x,a,b), add=T, col="red", lwd=2)
hist(Y, freq=F, main=paste("Beta from R:","a=",a,"b=",b), nclass=50, ylim=c(0,3))
curve(dbeta(x,a,b), add=T, col="red", lwd=2)