Wolfeyes Bioinformatics beta

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

練習問題 3.1 正規・コーシー-ベイズ推定量のモンテカルロ積分を2通りの方法で計算する

正規・コーシー-ベイズ推定量について， $x = 0,2,4$ だとしてモンテカルロ積分を計算する．

$\delta(x) = \frac{\int_{-\infty}^{\infty} \frac{\theta}{1+\theta^2}e^{-(x-\theta)^2/2} d\theta}{\int_{-\infty}^{\infty} \frac{1}{1+\theta^2}e^{-(x-\theta)^2/2} d\theta}$

さて，この式だけだと何のことか全然分からないのだが，式中で消えている係数をしっかり書くと，

分子： $\int_{-\infty}^{\infty}\theta\times\frac{1}{\pi(1+\theta^2)}\frac{1}{\sqrt{2\pi}}e^{-(x-\theta)^2/2}d\theta$
分母： $\int_{-\infty}^{\infty}\frac{1}{\pi(1+\theta^2)}\frac{1}{\sqrt{2\pi}}e^{-(x-\theta)^2/2}d\theta$

に分解され，結局はコーシー分布と正規分布の掛け算をしていることになる．模式的に書くと以下のようになる．

$\delta(x) = \frac{\int_{-\infty}^{\infty} \theta \times Cauchy \times Normal}{\int_{-\infty}^{\infty} Cauchy \times Normal}$

このように，上と下でそれぞれ積分しているので，モンテカルロ積分の実装においても上下別々に推定値を計算した後に $\delta(x)$ の推定値を求める．また，コーシー分布と正規分布を掛けた形となっているため，どちらの分布からも乱数を生成してモンテカルロ積分ができる．

ちなみに，本問題は4章の例4.2（P.107）にも同様の推定量を扱った問題が出てくる．

a.1 被積分関数をプロットする

左から順に，

$\delta(x)$ の分子の分布
$\delta(x)$ の分母の分布

を， $x = 0,2,4$ についてプロットしたのが以下の図になる．

delta_num <- function(t){t/(1+t^2)*exp(-(x-t)^2/2)}
delta_den <- function(t){1/(1+t^2)*exp(-(x-t)^2/2)}

par(mfrow=c(3,2))
for(x in c(0,2,4)){
  curve(delta_num, from=-10, to=10, main=paste("numerator : x=",x))
  curve(delta_den, from=-10, to=10, main=paste("denominator : x=",x))
}

a.2 コーシー・シミュレーションにもとづくモンテカルロ積分を計算する

コーシー分布から乱数を生成してモンテカルロ積分をシミュレーションする．イテレーションごとの推定値と標準分散のプロットは例3.3のプロットと同じで，黒の線が推定値の推移，上下の黄色の線が推定値の誤差の範囲となっている．

モンテカルロ積分の結果としては

x=0のとき：-0.01243
x=2のとき：1.282
x=4のとき：3.389

となった．

Nsim <- 10^4
par(mfrow=c(1,3))
for(x in c(0,2,4)){
  theta <- rcauchy(Nsim)
  num <- theta * dnorm(theta, mean=x)
  den <- dnorm(theta, mean=x)

  n <- num[den != 0]
  d <- den[den != 0] # ゼロ除算を起こさないようにしている
  y <- n / d
  estint <- (cumsum(n)/(1:length(n))) / (cumsum(d)/(1:length(d)))
  esterr <- sqrt(cumsum((y - estint)^2))/(1:length(y))
  plot(estint, xlab="Iterations", main=paste("x = ",x,", Estimate= ",mean(estint)),
       type="l", lwd=2, xlim=c(0,10^4), ylim=c(x-1,x+1))
  lines(estint+2*esterr, col="gold", lwd=2)
  lines(estint-2*esterr, col="gold", lwd=2)
}

b. 収束する様子を推定値の標準誤差でモニタリングする

上図の黄色の幅で示してある．

c. 正規近似によるモンテカルロ積分を計算する

以下のグラフでは，コーシー乱数のシミュレーションのグラフを重ねてプロットすることで，コーシー乱数を使った場合と正規乱数を使った場合の標準誤差を比較をしている．正規乱数を使ったシミュレーションの結果は，赤とピンクの線で表している．このプロットを見ると，正規乱数を使った場合の方がより誤差が少なく推定値を求められることがわかる．

黒＆黄色：コーシー乱数（問題 a.）
赤＆ピンク：正規乱数（問題 c.）

Nsim <- 10^4
par(mfrow=c(1,3))
for(x in c(0,2,4)){
  theta <- rnorm(Nsim, mean=x)
  num <- theta * dcauchy(theta)
  den <- dcauchy(theta)

  n <- num[den != 0]
  d <- den[den != 0]
  y <- n / d
  estint <- (cumsum(n)/(1:length(n))) / (cumsum(d)/(1:length(d)))
  esterr <- sqrt(cumsum((y - estint)^2))/(1:length(y))
  plot(estint, xlab="Iterations", main=paste("x = ",x,", Estimate= ",mean(estint)),
       type="l", lwd=2, col="red", xlim=c(0,10^4), ylim=c(x-1,x+1))
  lines(estint+2*esterr, col="pink", lwd=2)
  lines(estint-2*esterr, col="pink", lwd=2)
}

練習問題 3.9

この問題は，練習問題3.1の続き．

a.コーシー候補にもとづく受理・棄却アルゴリズムから推定量を計算する

$X \sim \mathcal{N}(\theta,1)$ $\theta \sim \mathcal{C}(0,1)$ において，事後分布 $p(\theta|x)$ は以下のように表される．

$p(\theta|x)=\frac{p(x|\theta)p(\theta)}{p(x)}$

ここで， $p(x|\theta) = \frac{1}{\sqrt{2\pi}}e^{-\frac{(x-\theta)^2}{2}}$ より，

$\begin{align} p(\theta|x) =& \frac{p(\theta)}{p(x)}\frac{1}{\sqrt{2\pi}}e^{-\frac{(x-\theta)^2}{2}} \ \propto& p(\theta)e^{-\frac{(x-\theta)^2}{2}} \end{align}$

となる．受理・棄却アルゴリズムにおける目標分布と候補分布は

目標分布 $f(\theta)$ ： $p(\theta)e^{-(x-\theta)^2/2}$
候補分布 $g(\theta)$ ： $\frac{1}{\pi(1+\theta^2)}$
受理・棄却条件： $U\times M \leq \frac{f(Y)}{g(Y)}$ ，（ $Y \sim \mathcal{C}(0,1)$ ， $U \sim \mathcal{U}_{[0,1]}$ ）

のようになる．受理・棄却アルゴリズムの実験は，練習問題3.1と同様に $x = 0,2,4$ で実行している．以下の図は，受理された乱数のヒストグラムと，その乱数の平均値の値を赤い線で示してある．

Nsim <- 10^4
x <- rep(0, Nsim)
M <- pi
f <- function(theta){ exp(-(xx-theta)^2/2)/(1+theta^2)}

par(mfrow=c(3,1))
for(xx in c(0,2,4)){
  for(i in 1:Nsim){
    y <- rcauchy(1)
    while(runif(1)*M > f(y)/dcauchy(y)){
      y <- rcauchy(1)
    }
    x[i] <- y
  }
  hist(x, freq=F, nclass=150, main=paste("x = ",xx))
  abline(v=mean(x), col="red", lwd=2)
}

b. 正規・コーシー-ベイズ推定量のモンテカルロ積分における，分母と分子の乱数の選択

分母と分子で同じ乱数を使う場合と違う乱数を使う場合で，どのような誤差の影響が出るかを比較する．今回は正規乱数を用いたモンテカルロ積分で試してみる．

黒＆黄色：分母と分子で同じ乱数 $\theta_i$ を使う場合
青＆水色：分母と分子で違う乱数 $\theta_i$ を使う場合

以下の図を見ると，xの値が大きい時に違う乱数を使う場合の方が標準誤差が大きい．つまり，分母と分子で別々にシミュレーションを行う際にも，乱数から生成する値を統一することで推定値の分散を押さえることができる．逆に言えば，分母と分子で別々の乱数を使うと，分母と分子どちらかが極端な値を引いてしまったときには推定値が変な値を取ってしまう場合がある．

Nsim <- 10^4
par(mfrow=c(1,3))
for(x in c(0,2,4)){
  theta <- rnorm(Nsim, mean=x)
  theta_2 <- rnorm(Nsim, mean=x)

  num <- theta * dcauchy(theta)
  den <- dcauchy(theta)
  num_2 <- theta * dcauchy(theta)
  den_2 <- dcauchy(theta_2)

  n <- num[den != 0]
  d <- den[den != 0]
  y <- n / d
  n_2 <- num_2[den_2 != 0]
  d_2 <- den_2[den_2 != 0]
  y_2 <- n_2 / d_2

  estint <- (cumsum(n)/(1:length(n))) / (cumsum(d)/(1:length(d)))
  esterr <- sqrt(cumsum((y - estint)^2))/(1:length(y))
  estint_2 <- (cumsum(n_2)/(1:length(n_2))) / (cumsum(d_2)/(1:length(d_2)))
  esterr_2 <- sqrt(cumsum((y_2 - estint_2)^2))/(1:length(y_2))

  plot(estint, xlab="Iterations", ylab="" ,type="l", lwd=2, col="black", xlim=c(0,10^4), ylim=c(x-1,x+1))
  lines(estint+2*esterr, col="gold", lwd=2)
  lines(estint-2*esterr, col="gold", lwd=2)

  par(new=T)
  plot(estint_2, xlab="Iterations", ylab="", type="l", lwd=2, col="blue", xlim=c(0,10^4), ylim=c(x-1,x+1))
  lines(estint_2+2*esterr_2, col="skyblue", lwd=2)
  lines(estint_2-2*esterr_2, col="skyblue", lwd=2)
}

参考

http://waxworksmath.com/Authors/N_Z/Robert_N_Casella/robert_n_casella.html

1章から3章途中までの解説とRのコードがある．今回の作図に関しては，このページのBook Notesを参考にさせていただいた．

Kindleに至るまでの小売時代のアマゾンを総括する1冊

「Aamazonは一体何の会社なのか」という問いに，ある人は本屋と答えるだろうし，またある人は何でも買えるネットショッピングサイトだと答えるだろう．インターネットの技術屋に聞けばAmazonはクラウドやAWSだと答えるだろうし，最近ではAmazonといえばKindleだと答える人が多いだろう．ジェフ・ベゾスが起こしたAmazonという会社は，今やインターネットでのショッピングのサービスを提供するという枠を乗り越え，AppleやMicrosoftなどに続いて不動の地位とブランドを獲得している．本書は，そのようなAmazonとジェフ・ベゾスが辿ってきた歴史を紐解いて，いかにAmazonが成長しジェフ・ベゾスが成功に導いてきたかに関する本である．

本書はジェフ・ベゾスの人生を軸に書かれている．まずはじめに，Amazonが成功するに至った2つの要因である

顧客を第一に考えること
長期的な視点をもって決断をすることに

の2つにまつわる話が描かれ，彼の作り上げた作品が紹介される．時には消費者の利便性を最大限に引き出すようなサービスを提供し，時には非情なまでに自動化されたシステムを作り上げる．輝かしい成功の影で行われている生々しい特許論争に関しても取り上げられる．そのようなAmazonをすべて統括してきたジェフ・ベゾスという人物が，いかにアントレプレナーとして優れており，その才能を発揮してきたかということが，彼の幼少時代から順に語られていく．彼はすぐにアントレプレナーとして頭角を現したわけでは無かったが，彼が大学を出た後に入った会社では非情に優れた技術者として出世していったという逸話や，インターネットを使った書籍販売というアイデアを思いつくまでの経緯などからも，Amazonに繋がる一連の流れをうかがい知ることができる．Amazonの会社を起こしてからは，今からは想像も出来無いようなベンチャーとしての苦労話や逸話が出てくる．そして最後の数章では，現在のAmazonのイノベーションの中心とも言えるKindleの誕生とその戦略にフォーカスが移る．

個人的な感想としては，興味はやはりAmazonが現在推し進めているKindleやAWSに集中してしまう．ネットショッピングとしてのAmazonはもはやネットに当然のようにあるものとして，空気のような必要不可欠なものであり，普段はその重要性や革新性に意識を向けないような存在になってしまった．それはまさに，AppleがiPodを開発して音楽業界に進出したときのような時代の分かれ目を感じさせる．そういった現状を踏まえると，本書は小売時代のAmazonという企業の歴史を総括する1冊として上手く収まっているといえる．その反面，本書におけるKindleの扱いはどうしても付け加え感が拭えず，本書は全体的に伝記的な作品とも言いにくくマーケティングや経営戦略的な部分も薄いといった難点はあるものの，これからの電子書籍戦争を見守る中でAmazonの立ち位置が分かりやすく示されているところはとても評価できる部分だ．本書で見られるジェフ・ベゾスの一貫した理念を考えると，書籍から電子書籍へと媒体が変わったとしても恐らく根本の部分では何も変わらないのだろう．

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

内容：乱数を利用して積分を数値的に求める

モンテカルロ積分とは，乱数を利用して単変量や多変数の積分を近似し，解析的に解けない問題を数値的に解くという手法である．これは2章で示したように，任意の分布に従う乱数を理論上無限に生成できるからこそ成り立つ手法であり，本手法はi.i.dを前提とした大数の法則や中心極限定理のような確率論に落としこんで解析している．具体的なモンテカルロ法の利用例としては確率分布を積分といったことが挙げられ，ベイズなどに限らず確率分布を仮定するような統計手法などにおいて必要となる．

例 3.3 モンテカルロ積分の具体例

図の上半分は $h(x) = [\cos(50x) + \sin(20x) ]^2$ を0から1の間で図示したもの．モンテカルロ積分ではこの曲線の下側の面積を求めることになる．図の下半分は，実際にモンテカルロ積分のシミュレーションをしたもの．横軸がシミュレーション回数（サンプル数m）で，縦軸が経験平均 $\bar{h}_m = \frac{1}{m}\sum_{j=1}^{m}h(x_j)$ を示したもの．黒の線が平均値の推移，上下の黄色の線が推定標準誤差に基づく信頼幅，赤の直線がRのintegrate関数で求めた積分値を示している．この図から，シミュレーション回数を重ねるごとに数値的に求めた値が実際の積分値に収束していく様子がわかる．

h <- function(x){(cos(50*x) + sin(20*x))^2}
par(mar=c(2,2,1,1), mfrow=c(2,1))
curve(h, xlab="Function", ylab="", lwd=2)
integrate(h,0,1)
x <- h(runif(10^4))
estint <- cumsum(x)/(1:10^4)
esterr <- sqrt(cumsum((x-estint)^2))/(1:10^4)
plot(estint, xlab="Mean and error range", type="l", lwd=2, ylim=mean(x)+20*c(-esterr[10^4], esterr[10^4]))
lines(estint+2*esterr, col="gold", lwd=2)
lines(estint-2*esterr, col="gold", lwd=2)
abline(h=integrate(h,0,1)$value,lwd=2,col="red")

例 3.4 正規分布のモンテカルロ積分における値の精度と効率

「裾」と呼ばれる部分は正規分布において釣鐘型の山の頂点から遠く離れた部分のことを指しており，正規分布の平均から外れるということは起こる確率がきわめて低い部分である．そのため，正規乱数を基にしたモンテカルロ積分においては，山の頂点に近い部分の値は乱数でよく引きあてるので数は十分なのだが，裾の部分の「値が大きい/小さい」値は乱数でなかなか引き当てることができない．精度に大きな影響を与える部分でありながらもシミュレーションを重ねないと十分な数を得ることができないため，古典的なモンテカルロ積分の方法で精度を上げるためには，試行回数を多くしなければならない．

以下の実験では，実際に有効数字4ケタで正確な数値を出すために，10⁸ 回ものシミュレーションを行なっている．最終的な表3.1は，行が試行回数，列が $x \sim \mathcal{N}(0,1)$ のときのxに対応している．一番下の行の値がboundという変数のそれぞれ値と同じになっていることから，10⁸ もの試行回数が必要であることがわかる．

x <- rnorm(10^8)
bound <- qnorm(c(0.5, 0.75, 0.8, 0.9, 0.95, 0.99, 0.999, 0.9999))
res <- matrix(0, ncol=8, nrow=7)
for(i in 2:8)
  for(j in 1:8)
    res[i-1, j] <- mean(x[1:10^i] < bound[j])
matrix(as.numeric(format(res, digi=4)), ncol=8)

# > matrix(as.numeric(format(res, digi=4)), ncol=8)
# [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]
# [1,] 0.5200 0.7600 0.8000 0.9300 0.9800 1.0000 1.0000 1.0000
# [2,] 0.4760 0.7400 0.7880 0.8980 0.9550 0.9950 1.0000 1.0000
# [3,] 0.4980 0.7477 0.7967 0.8981 0.9480 0.9889 0.9993 0.9999
# [4,] 0.4994 0.7503 0.7996 0.8992 0.9504 0.9899 0.9989 0.9999
# [5,] 0.4999 0.7503 0.8002 0.9001 0.9501 0.9899 0.9990 0.9999
# [6,] 0.4999 0.7500 0.8000 0.9001 0.9500 0.9900 0.9990 0.9999
# [7,] 0.5000 0.7500 0.8000 0.9000 0.9500 0.9900 0.9990 0.9999

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

練習問題 2.18

a. f(x)とMg(x)をプロットする

定数Mは例2.7のようにoptimize()関数で求めている．以下の図はfとgの分布を図示したもので，黒い曲線fに対して，青い曲線Mg(x)が全体を覆い被さるようになっていることがわかる．

f <- function(x){exp(-x^2/2) * (sin(6*x)^2 + 3*cos(x)^2*sin(4*x)^2 + 1)}
g <- function(x){exp(-x^2/2)/sqrt(2*pi)}
M <- optimize(f=function(x){f(x)/g(x)}, interval=c(0,1),maximum=T)$objective

par(mfrow=c(1,1))
curve(f(x),from=-5,to=5,ylim=c(0,5))
curve(M*g(x),add=T,col="blue")

b. 受理・棄却アルゴリズムを使ってfから2500個の乱数を生成する

ここで標準正規分布からの乱数を使うのは，g(x)が標準正規密度だから（P.59の受理・棄却法のアルゴリズムの1.にある $Y \sim g$ の部分に該当する）．

この問題では2500個の受理された乱数が欲しいので，P.61のコードを真似て作成している．ただし，このコードは少々効率が悪く，決められた数の乱数を作ったらループが終わりというものではない（受理確率を元にして多めに作ってる）ので，ヒストグラムとして図示する際には2500個のみを取り出して作図している．

Nsim <- 2500
x <- NULL
while(length(x) < Nsim){
  y <- rnorm(Nsim*M)
  x <- c(x, y[runif(Nsim*M,0,M) * g(y) < f(y) ])
}
hist(x[1:Nsim],freq=F,nclass=150,main="Nsim = 2500")
curve(f(x)/integrate(f,-Inf,Inf)$value, add=T, col="red")
# curveの中のintegrateに関しては後述(本問のc.)

概ねfの分布に沿った乱数が生成されているようだ．乱数が2500個だけでは少し不明瞭のため，Nsimの数を10⁴ にして再試したのが下の図．

c. シミュレーションから得られた受理率から正規化定数を求める

P.60にあるように，受理確率は基本的に1/Mだけれども，正規化されていない関数に関しては定数CがMに吸収されているので注意が必要となる．今回の場合，シミュレーションから求めた受理確率rを使うことで，以下の式から正規化定数を近似することができる．

$r = \frac{1}{C \times M}$

今回の場合，Mは10.94，rは0.54となったので，Cは0.17程度だと見積もった．

Nsim <- 10^4
x <- NULL # xは受理した乱数
z <- NULL # zは生成した乱数すべて
while(length(x) < Nsim){
  y <- rnorm(Nsim*M)
  x <- c(x, y[runif(Nsim*M,0,M) * g(y) < f(y) ])
  z <- c(z, y)
}
r <- length(x)/length(z)
C <- 1/(M*r)

> M; r; C
[1] 10.94031
[1] 0.5369048
[1] 0.1702445

さて，答え合わせ(?)だが，正規化定数は単純にfの $[-\infty,\infty]$ の積分を求めることによって計算することができる．これにはRのintegrate()関数を使って-infからinfまでを計算すればよい．確率分布の定義から，これが1になるように正規化定数を定めれば良いということで，以下のように計算した結果，0.1696543となった．

> 1/integrate(f,-Inf,Inf)$value
[1] 0.1696543

練習問題 2.19 二重指数分布が候補分布の受理・棄却アルゴリズムにおいて標準正規分布から乱数を生成する際のMとg(x)の最適化

$f(x) = \frac{1}{\sqrt{2\pi}} \exp(-\frac{x^2}{2}) \ g(x) = \frac{\alpha}{2} \exp(-\alpha |x|)$

$f/g = \frac{\sqrt{2/\pi}}{\alpha} \exp(\alpha|x|-\frac{x^2}{2})$ より，上界Mを取るときのx_maxの値を調べるには，とりあえずxで微分して0になるときの値を求めれば良い．expの部分を微分すると $(\alpha - x)$ みたいな項が出てくるので， $x_{max} = \alpha$ （厳密には $x_{max} = \pm \alpha$ ）．

次に，受理率を最適化するということは1/Mの値を大きくすれば良いので，結局はMの値の最大値を求めれば良い． $x_{max} = \alpha$ のときのMを $\alpha$ に関して微分して0になるときの値ということで計算していくと $\frac{(\alpha^2-1)\exp(\alpha^2/2)}{\alpha^2}$ みたいな項が出てくるので，つまり $\alpha = 1$ の時にMが最大値を取り，受理率が最適化される．

練習問題 2.22 切断正規分布から正規乱数を生成する

a. 切断正規分布から乱数を生成する

library(msm) # dtnorm()のため
mu <- 0
sigma <- 1

par(mfrow=c(4,2))
for(a in c(-2.32, -0.84, -0.67, 0, 0.67, 0.84, 2.32)){
  N <- 10^4
  X <- rep(0,N)
  for(i in 1:N){
    z <- rnorm(1, mean=mu, sd=sigma)
    while(z < a){
    z <- rnorm(1, mean=mu, sd=sigma)
    }
  X[i] <- z
  }
  hist(X, freq=F, nclass=150, main=paste("Truncated Normal Distribution: a = ",a))
  curve(dtnorm(x, mean=mu, sd=sigma, lower=a, upper=Inf), add=T, col="red", lwd=2)
}

b. シミュレーションので受理する確率を求める

上で示したアルゴリズムで考える．z < aならアタリ，それ以外ならハズレというようなときにどれだけ数を撃てばz<aに入るかということを考えると，正規分布のz < aの面積を求めればいいことがわかる．なお，ここで出てくる $\Phi(x)$ という関数は定義されていないが，以下のような累積標準正規分布（本書P.78に出てくる）のことだと思われる．

$\Phi(t) = \int_{-\infty}^{t} \frac{1}{\sqrt{2\pi}}e^{-y^2/2}dy$

この形式になるように，aからの累積正規分布 $F(x|\mu, \sigma^2, a) \sim \int_{a}^{\infty} \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}dx$ を変形していけばよい．

ここで， $v = -\frac{x-\mu}{\sigma}$ とおくと， $dv = -\frac{dx}{\sigma}$ ， $x:a \rightarrow \infty$ のとき， $v:-\infty \rightarrow \frac{\mu-a}{\sigma}$ をとる．よって，

$F(v|\mu,\sigma^2,a) = \int_{\infty}^{\frac{\mu-a}{\sigma}} \frac{1}{\sqrt{2\pi}}e^{-\frac{v^2}{2}}dv = \Phi(\frac{\mu-a}{\sigma})$   となる．

また，aの値が裾にある（値が大きくなる）と， $\Phi(\frac{\mu-a}{\sigma})$ の値が小さくなるので，シミュレーション回数が多くなる．

c. 切断正規分布の乱数を正規分布から生成することを考える

目標分布； $\exp(-(x-\mu)/2\sigma^2)\mathbb{I}_{x \geq a}$
候補分布： $\frac{1}{\sqrt{2\pi}}\exp(-\frac{(x-\bar{\mu})^2}{2})$

より

$\frac{f(x)}{g(x)}=\frac{\exp(-(x-\mu)/2\sigma^2)\mathbb{I}_{x\geq z}}{\frac{1}{\sqrt{2\pi}}\exp(-\frac{(x-\bar{\mu})^2}{2})}\$ $= \sqrt{2\pi}\exp(-x\bar{\mu}+\frac{1}{2}\bar{\mu}^2)\mathbb{I}_{x \geq a}$

となり候補分布の制約を満たすので，正規候補にもとづく受理・棄却アルゴリズムは切断正規分布から乱数を生成するために使える．

次に，このときの最適なMを考える．x=aにおいて $\frac{f(a)}{g(a)} = \sqrt{2\pi}\exp(-a\bar{\mu}+\frac{1}{2}\bar{\mu}^2)\mathbb{I}_{x \geq a}$ より，この式を $\bar{\mu}$ で微分すると

$\Bigl(\frac{f(a)}{g(a)}\Bigr)' = \sqrt{2\pi}(-a+\bar{\mu})\exp(-a\bar{\mu}+\frac{1}{2}\bar{\mu}^2)\mathbb{I}_{x \geq a}$

となるので，Mは $\bar{\mu}=a$ のとき最適値をとる．

d. 候補分布に指数分布を使うことを考える

目標分布； $\exp(-(x-\mu)/2\sigma^2)\mathbb{I}_{x \geq a}$
候補分布： $g_{\alpha}(z)=\alpha e^{-\alpha(z-a)}\mathbb{I}_{x \geq a}$

$\begin{align}\frac{f(x)}{g(x)}&=\frac{\exp(-(x-\mu)/2\sigma^2)\mathbb{I}_{x\geq a}}{\alpha e^{-\alpha(z-a)}}\&=\frac{1}{\alpha}\exp(-\frac{x^2}{2}+\alpha(x-a))\mathbb{I}_{x \geq a}\end{align}$

より，

$\frac{f(x)}{g(x)} \propto e^{\alpha(x-a)}e^{-z^2/2}$

となる．上式を前と同じように微分すると $(-x+\alpha)$ が出てくるので， $x = \alpha$ のとき最大値を取り， $\exp(\alpha^2/2-\alpha a)$ で極限が与えられる．

次に， $a = \alpha$ の場合に正当な候補密度になることを導く．f/gの比が0より大きいことを示せば良いのだが，上式を見ると結局expの外側の $\frac{1}{\alpha}$ ，つまり $\alpha$ を調べないといけない．ということで， $\alpha$ で微分して最適値の時の値を求める．

$x = \alpha$ のとき $\frac{f(\alpha)}{g(\alpha)} = \frac{1}{\alpha}\exp(\frac{\alpha^2}{2}-\alpha a)$ なので，

$\Bigl(\frac{f(\alpha)}{g(\alpha)}\Bigr)' = -\frac{1}{\alpha^2}\exp(\frac{\alpha^2}{2}-\alpha a)+\frac{1}{\alpha}(\alpha-a)\exp(\frac{\alpha^2}{2}-\alpha a)$ $=\frac{1}{\alpha}\bigl(-\frac{1}{\alpha^2}+(\alpha-a)\bigr)\exp(\frac{\alpha^2}{2}-\alpha a)$

より

$-\frac{1}{\alpha^2}+(\alpha-a) = 0$

$\alpha^2 - a \alpha - 1 = 0$

よって， $\alpha = \frac{a \pm \sqrt{a^2 + 4}}{2}$ となる．

じゃあプラスかマイナスかどっちなんだという話になるのだが，ここで $\alpha \geq a$ という前提条件を使うと，

$\alpha = \frac{a \pm \sqrt{a^2 + 4}}{2} \geq a$

が成り立つ． $\sqrt{a^2 + 4}$ というのはaよりちょっとだけ大きい数なので，上式を満たすaは $\alpha = \frac{a + \sqrt{a^2 + 4}}{2}$ となり，この式は常に正である．よって， $a = \alpha$ の場合に，正当な候補密度になる．

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

内容：逆変換法や一般変換法で生成できない分布から乱数を生成する

逆変換法や一般変換法では，一様分布から生成される乱数に何らかの関数を通すことで任意の乱数を生成してきた．しかし，このような直接的な方法で乱数を生成できない分布の場合には，間接的な方法によって乱数を生成する必要がある．そのような場合には，本来求めたい乱数の分布とは違った分布から乱数を生成し，それが受理できる値か棄却できる値かを選り分けることで，間接的に求めたい分布から生成された乱数を表現する．

目標密度（target）: $f$
手段密度・候補密度（candidate）: $g$

この受理・棄却法では1つの乱数の値 $X$ につき，2つの乱数 $Y \sim g$ と $U \sim \mathcal{U}_{[0,1]}$ を考える．前者は候補密度から乱数を生成したもの，後者は受理・棄却に用いる一様乱数．この2つの乱数を，以下の式にあてはめて受理するか棄却するかを判断する．

$U \leq \frac{f(Y)}{M \times g(Y)}$

ここで登場するMは提案分布と候補密度から見積もった定数で，候補密度の制約の中で登場する値．基本的にある値より大きければ何でも良いのだが，シミュレーションの効率に関わってくる．これが小さい値であればあるほど，シミュレーションで棄却される乱数の候補が少なくて済むので，結果として効率が良くなる．最適なMは以下の式で与えられる（UserR!資料pdf p.51）．

$M = \sup_x \frac{f(x)}{g(x)}$

例 2.7, 例2.8 ベータ分布の乱数を生成する(2通り)

黒の点が一様分布/ベータ分布から生成した乱数，赤の点がアルゴリズムで受理された乱数（今回求めたかった乱数），青の線が極限 $M g(x)$ ．

受理・棄却の条件

例2.7: ，，，
- ここでgは何らかの分布ではなく「gは1に等しい」としているので，上式の1つ目の式の中にgが含まれていない．この場合，gで生成される乱数は一様乱数．
- 一様分布のグラフを考えてみるとわかりやすい．区間[0,1]の間で常に値1を取る関数なので，上では1に等しいと表している．
例2.8: $Ug(Y) < f(Y)$ ， $U \sim \mathcal{U}_{[0,M]}$ ， $f \sim \mathcal{Be(2.7,6.3)}$ ， $g \sim \mathcal{Be(2,6)}$

例2.8における「提案分布」という言葉は，新しく設定した候補分布のこと．

# 例2.7
Nsim <- 2500
a <- 2.7
b <- 6.3
M <- 2.67
u <- runif(Nsim, max=M)
y <- runif(Nsim)
x <- y[u < dbeta(y,a,b)]

y2 <- rbeta(Nsim,2,6)
u2 <- runif(Nsim, max=1.67) * dbeta(y2,2,6)
x2 <- y2[u2 < dbeta(y2,2.7,6.3)]

par(mfrow=c(1,2))
plot(y, u, xlim=c(0,1.0), ylim=c(0,2.8), cex=0.5, xlab="y", ylab="u.g(y)")
par(new=T)
plot(x, u[u < dbeta(y,a,b)], col="red", xlim=c(0,1.0), ylim=c(0,2.8), cex=0.5, xlab="y", ylab="u.g(y)")
curve(dbeta(x,a,b), add=T, col="red", lwd=2)
abline(h=M, lwd=2, col="blue")

plot(y2, u2, xlim=c(0,1.0), ylim=c(0,5), cex=0.5, xlab="y", ylab="u.g(y)")
par(new=T)
plot(x2, u2[u2 < dbeta(y2,2.7,6.3)], col="red", xlim=c(0,1.0), ylim=c(0,5), cex=0.5, xlab="y", ylab="u.g(y)")
curve(dbeta(x,2.7,6.3), add=T, col="red", lwd=2)
curve(dbeta(x,2,6)*1.67, add=T, col="blue", lwd=2)

上の図の解釈

乱数の値が赤の分布の線より中に入れば受理となる．候補として生成される乱数は青の分布の線の中でランダムに生成されるので，赤と青の線の間の部分が棄却された乱数になる．この赤と青の線の間の部分が少なければ少ないほど，この受理・棄却法の効率も良くなる．じゃあ効率を良くするにはどうすればいいかというと，いじることのできるポイントは2つ

提案分布gを良い感じにする（青の線の形を決める）
定数Mを出来るだけ小さくする（青の線の高さを決める，fとgの間の幅を決める）

この受理・棄却法の前提として目標分布（赤の線）から乱数は生成できないのだから，100％受理されるようなシミュレーションは理論上不可能．だから，乱数が生成できるような候補分布を使って，定数Mを定めて，間接的に乱数を発生させようということ．

練習問題 2.7 目標分布fと候補密度gがともにベータ分布だった場合に，どのような制約があるのか

$a \leq \alpha, b \leq \beta$ が必要なことの証明

$f(x) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} x^{\alpha-1} (1-x)^{\beta-1}$ ， $g(x) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} x^{a-1} (1-x)^{b-1}$ より

$\frac{f(x)}{g(x)} = \frac{\Gamma(\alpha+\beta)\Gamma(a)\Gamma(b)}{\Gamma(a+b)\Gamma(\alpha)\Gamma(\beta)} x^{\alpha-a} (1-x)^{\beta-b}$

となり，比 $f/g$ が有界になる（どっかのxの値が $\inf$ になったりしない）には，xや(1-x)のベキ数の部分が0以上でである必要がある．よって $a \leq \alpha, b \leq \beta$ が必要．

上の細かい証明は省略するが，要するにベータ分布 $\mathcal{Be}(a,b)$ において，aとbが1未満の時には0と1の両端で値が跳ね上がる．以下の図は，「パターン認識と機械学習上」のP.70を参考にベータ分布を作図したもの．これを見ると，aとbが1のときにちょうど一様分布のような形をとり，それより小さいと下に凸，それより大きいと上に凸のような分布の形になることがわかる．なので，ベータ分布が有界であるにはxや(1-x)のベキ数の部分がポイントになってくる．

# 練習問題 2.7 ベータ関数の分布を図示してみる
par(mfrow=c(2,2))
curve(dbeta(x,0.1,0.1),col="red",lwd=2,ylim=c(0,3),main="Be(0.1,0.1)")
curve(dbeta(x,1,1),col="red",lwd=2,ylim=c(0,3),main="Be(1,1)")
curve(dbeta(x,2,3),col="red",lwd=2,ylim=c(0,3),main="Be(2,3)")
curve(dbeta(x,8,4),col="red",lwd=2,ylim=c(0,3),main="Be(8,4)")

$a = \left\lfloor \alpha \right\rfloor$ ， $b = \left\lfloor \beta \right\rfloor$ の証明

$f(x)/g(x)$ を微分して最大値を取る時のxを求めると， $x_{max} = \frac{\alpha-a}{\alpha-a+\beta-b}$ となる．よって，上で証明した $a \leq \alpha, b \leq \beta$ を満たしつつ自然数aとbが最適な値を取るには， $a = \left\lfloor \alpha \right\rfloor , b = \left\lfloor \beta \right\rfloor$ となる．これもまあ厳密な証明はできなくても直感的に考えれば明らかだろう．

練習問題 2.8 二重指数分布から正規乱数を生成する

$f(x) = \frac{1}{\sqrt{2\pi}}\exp(-\frac{x^2}{2})$
$g(x|\alpha) = \frac{\alpha}{2}}\exp(-\alpha |x|)$

a. 定数Mの取りうる範囲と最小値を求める

$\frac{f(x)}{g(x|\alpha)} = \sqrt{\frac{2}{\pi}}\alpha^{-1}\exp(-\frac{x^2}{2} + \alpha|x|)$

$x = \alpha$ のとき， $\frac{f(\alpha)}{g(\alpha)} = \sqrt{\frac{2}{\pi}}\alpha^{-1}\exp(-\frac{\alpha^2}{2} + \alpha^2) = \sqrt{\frac{2}{\pi}}\alpha^{-1}e^{\alpha^2/2}$ となり，与式を満たす．

Mの最小値は，上式を微分すると $(-\alpha^2+1)\exp(-\frac{\alpha^2}{2})$ が出てくるので， $\alpha = \pm 1$ となる．候補分布に課された制約(P.59)より，fとgどちらか一方だけ負になることは無いので， $\alpha = 1$ が正解．（二階微分すれば-1のとき最小値，1のとき最大値だとわかる）