棄却法
密度πをもつ確率分布からの乱数が発生したいとします。
πに対して、そのサポート(π(x)>0となるようなxの範囲のこと)を包むようなサポートを持っている密度qをもつ分布からの乱数が発生できるとします。
ここで、πのサポート上の任意のxについて、以下のような定数Cが存在するとします。
このような定数Cについて、以下のアルゴリズムにより発生される乱数xは、密度πをもつ分布に従っています。
[アルゴリズム：棄却法]
①qから乱数yを発生させる
②(0,1)の一様乱数uを発生させる
③ であれば、このyを乱数xとする。ステップ①に戻る。
ちなみに、棄却法に限らず、このような我々が乱数を発生させることができる(簡単な)分布qのことを、提案分布と呼びます。乱数が欲しい分布πのことを目的分布と呼ぶこともあります。
なお、棄却法による、πに従う乱数を1個発生させるために必要な、qから発生させる乱数の個数は平均C個ですので、このような値Cは可能な限り低くとると、効率が良いです。
ちなみに、πに対してqの関数の減衰が明らかに速いと、π(x)/q(x)を上から抑える数Cはいくらでも大きくなってしまうので、注意。
理論的な背景は、鎌谷『モンテカルロ統計計算』の58p～や、小澄『ベイズ計算統計学』の30p～を読むと詳しいので、読んでみてください。
本当かよ？
さて、今までの説明は本当でしょうか？ 私が嘘をついている可能性がありますね。
ま、試しにちょっとやってみましょうか。
位置パラメータ0、尺度パラメータ1、形状パラメータ3のフレシェ分布を考えてみましょ～
curve(3*x^(-3-1)*exp(-x^(-3)) , xlim=c(0,7) , main = "フレシェ分布(形状パラメータ=3)")
ちなみにこいつはタイプII極値分布と呼ばれる分布になっていて、極値統計学との関係性が深いのですが、まあ今はそこら辺の話は置いときましょう。
この分布からの乱数を考えてみます。
こいつは結構裾が厚いので、例えば同じく正の値のみをとる分布として指数分布を考えてみると、
これについて、指数分布の減衰の方が圧倒的に早いので、π(x)/q(x)は、特に裾の方で以下のようになり、
curve(3*x^(-3-1)*exp(-x^(-3)) , xlim=c(5,15) , main = "フレシェ分布と指数分布" , col="red")
curve(exp(-x) , add=T , col="blue")
legend("topright" , legend = c("フレシェ分布" , "指数分布") , col=c("red" ,"blue") , lty=1,cex=1.5 , lwd=2)
curve(3*x^(-4)*exp(x-x^(-3)) , xlim=c(0,10) , main="比率")
これを上から抑えるような定数は存在しないので、棄却法は実行できませんね。
じゃあ、もう仕方ないんで、私たちが良く知っている最強に裾が厚いアイツを使いましょう。
(標準)コーシー分布を提案分布として利用します。
これとフレシェ分布の密度との比率は以下のようになります。
curve(3*pi*(x^(-2)+x^(-4))*exp(-x^(-3)) , xlim=c(0,10) , main="比率")
というわけで、コーシー分布を提案分布として、棄却法アルゴリズムを実行してみます。
図から見る限り、まあ大体7くらいで抑えられるみたいなんで、C=7.2として実行してみましょう。
res <- c(0)
N <- 10000
flag <- 0
for(i in 1:N){
y <- rcauchy(n=1)
u <- runif(1)
if(u < (3*y^(-4)*exp(-y^(-3)))/(7.2*1/(pi*(y^2+1))) && y>0){
res <- append(res , y)
}else{
flag <- flag+1
}
}
res
flag
さて、手に入れた乱数resをプロットしてみましょう。
plot(density(res) , xlim = c())
curve(3*x^(-3-1)*exp(-x^(-3)) , col="red" , add=T)
legend("topright" , legend="フレシェ分布" , lty=1 , lwd=2 , col="red" , cex=1.7)
このように、棄却法により手に入れた乱数が、フレシェ分布に従っていそうだと確認できました。