- n人の投擲選手が居て、その人たちの投擲記録について賭けをすることを考える
- n人の投擲選手の記録はガンマ分布に従うと予想する。ガンマ分布を選ぶのは、0以上の１ピークの分布であって、まあまあいい感じのモデルだから
- ある選手の記録は以下のような確率分布になる
x <- seq(from=0,to=100,length=101)
m1 <- 30
sc1 <- 5
sh1 <- m1/sc1
y1 <- dgamma(x,shape=sh1,scale=sc1)
plot(x,y1,type="l")
- 選手人数n=4として、それらの誰が投げるかは等確率で決まるものとする
- このとき、投擲記録の確率密度分布は以下のようになる
n <- 4
m <- c(20,40,50,60)
sc <- c(1,2,1,3)
sh <- m/sc
ys <- matrix(0,n,length(x))
for(i in 1:n){
ys[i,] <- dgamma(x,shape=sh[i],scale=sc[i])
}
par(mfcol=c(1,2))
matplot(x,t(ys),type="l")
y <- apply(ys,2,mean)
plot(x,y,type="l")
par(mfcol=c(1,1))
- さて、賭けというのは、こういう賭けである
- 賭ける人は、平均Mで標準偏差Sの正規分布を指定して良いことになっている。４選手の誰かが投げて出た記録に対して、指定した正規分布の確率密度を得点として得られる、というのが、この賭けのルールである
- 今から賭けようと思うが、正規分布のパラメタの組(M,S)について、得点の期待値が知りたい。さて、どうするか
- 地道にやってみよう
- xの値を0から小刻みに増やして、それぞれの生起確率を出して、そのときの得点を積み上げていくことにする
- それぞれのチェックポイントには生起確率と得点がある
M <- 40
S <- 20
n.r <- 10^4
x.check <- seq(from=0,to=1000,length=n.r)
Point <- 0
sum.prob <- 0
for(i in 1:length(x.check)){
tmp.prob <- 0
for(j in 1:n){
tmp.prob <- tmp.prob + dgamma(x.check[i],shape=sh[j],scale=sc[j])
}
tmp.prob <- tmp.prob/n
Point <- Point + dnorm(x.check[i],M,S) * tmp.prob
sum.prob <- sum.prob + tmp.prob
}
Point <- Point/sum.prob
Point
> Point
[1] 0.01479605
-
- これは、言い方を変えると、ある指定範囲に一様乱数のようなもの(等間隔指定点)を発生させて、その生起確率で重みをつけて、重み付き平均を出している
- 重み付き平均を期待値とするのはインポータンスサンプリングのやり方だから、「一様分布」と「確率密度分布」とで重みを決めて、乱点のようなものとして、等間隔指定点セットを取っている場合に相当する
- 生起確率密度分布から乱数を発生させて推定してみる
rs <- matrix(0,n,n.r/n)
for(i in 1:n){
rs[i,] <- rgamma(n.r/n,shape=sh[i],scale=sc[i])
}
hist(rs)
Point <- 0
for(i in 1:length(rs)){
Point <- Point + dnorm(rs[i],M,S) * 1/length(rs)
}
Point
> Point
[1] 0.01485484
-
- これは、まったく重みを付けずに、超地道に乱数を発生させて期待値を出している
- さて。インポータンスサンプリングでは、先に一様乱数のようなものを使ったが、もっと別の確率分布を使ってみよう、というもの
- どうしてそうするとよいか、というと、重み付き平均を出すときの、「重み」ｘ「値」が大きめな値で揃っているときに、期待値の推定精度が上がることが知られているから
- なので、何かうまい分布を取りたいな、という話
- 例えば、投擲距離の理論分布の平均と標準偏差を、上で生成した乱数のそれを取って、その平均と標準偏差を持つような正規乱数を発生させ、その正規分布とガンマ分布の混合分布である理論分布との比を重みづけにしてやってみる
mm <- mean(rs)
ss <- sd(rs)
rs <- rnorm(n.r,mm,ss)
Point <- 0
for(i in 1:length(rs)){
tmp.prob <- 0
for(j in 1:n){
tmp.prob <- tmp.prob + dgamma(rs[i],shape=sh[j],scale=sc[j])
}
tmp.prob <- tmp.prob/n
tmp.point <- dnorm(rs[i],M,S) * 1/length(rs)
tmp.weight <- tmp.prob/dnorm(rs[i],mm,ss)
Point <- Point + tmp.point * tmp.weight
}
Point
> Point
[1] 0.01480977
- ここまでで、変な確率密度分布をする確率変数があって、それに応じて、何かしらの得点を想定し、その期待値を求める話をしてきた
- 区間からのの等間隔点を出すのもやったが、その代わりに一様乱数を使うことにして、同様のことをやってみる
- ただし、乱点の数を少な目にして、モンテカルロ推定値がばらつくようにして、そのバラツキが大きいか小さいかを見てみる
M <- 40
S <- 20
n.r <- 100
n.iter <- 1000
Points.unif <- Points.gamma <- Points.norm <- rep(0,n.iter)
for(ii in 1:n.iter){
rs <- runif(n.r,0,1000)
Point <- 0
sum.prob <- 0
for(i in 1:length(rs)){
tmp.prob <- 0
for(j in 1:n){
tmp.prob <- tmp.prob + dgamma(rs[i],shape=sh[j],scale=sc[j])
}
tmp.prob <- tmp.prob/n
tmp.point <- dnorm(rs[i],M,S) * 1/length(rs)
tmp.weight <- tmp.prob/dunif(rs[i],0,1000)
Point <- Point + tmp.point * tmp.weight
}
Points.unif[ii] <- Point
rs <- matrix(0,n,n.r/n)
for(i in 1:n){
rs[i,] <- rgamma(n.r/n,shape=sh[i],scale=sc[i])
}
Point <- 0
for(i in 1:length(rs)){
Point <- Point + dnorm(rs[i],M,S) * 1/length(rs)
}
Points.gamma[ii] <- Point
rs <- rnorm(n.r,mm,ss)
Point <- 0
for(i in 1:length(rs)){
tmp.prob <- 0
for(j in 1:n){
tmp.prob <- tmp.prob + dgamma(rs[i],shape=sh[j],scale=sc[j])
}
tmp.prob <- tmp.prob/n
tmp.point <- dnorm(rs[i],M,S) * 1/length(rs)
tmp.weight <- tmp.prob/dnorm(rs[i],mm,ss)
Point <- Point + tmp.point * tmp.weight
}
Points.norm[ii] <- Point
}
Points.three <- cbind(Points.unif,Points.gamma,Points.norm)
boxplot(Points.three)
> apply(Points.three,2,mean)
Points.unif Points.gamma Points.norm
0.01498932 0.01479045 0.01477595
> apply(Points.three,2,sd)
Points.unif Points.gamma Points.norm
0.0068368258 0.0003667666 0.0006546649
-
- わかることは、生起確率密度そのものから乱数が発生できる場合(中央)はばらつきが小さく、一様乱数の場合はばらつきが大きいこと、正規乱数は結構、ばらつきが小さいことである
- このことから、一様乱数より正規乱数の方がよいことである。それは、正規乱数が一様乱数よりも、本当の生起確率分布と得点値との積の分布に近いから
- では、生起確率密度が高いときに得点が小さいような場合はどうなるだろうか
- このときに正規乱数の方を、「賭ける得点分布」にしてやってみる。なぜなら、正規乱数の密度ｘ得点に近いときが良いのだから、『得点が高い場合には、投擲記録の出る確率に大差がないような場合』には、得点の分布のみの高低が影響しそうだから
M <- 150
S <- 30
mm <- M
ss <- S
n.iter <- 1000
Points.unif <- Points.gamma <- Points.norm <- rep(0,n.iter)
for(ii in 1:n.iter){
rs <- runif(n.r,0,1000)
Point <- 0
sum.prob <- 0
for(i in 1:length(rs)){
tmp.prob <- 0
for(j in 1:n){
tmp.prob <- tmp.prob + dgamma(rs[i],shape=sh[j],scale=sc[j])
}
tmp.prob <- tmp.prob/n
tmp.point <- dnorm(rs[i],M,S) * 1/length(rs)
tmp.weight <- tmp.prob/dunif(rs[i],0,1000)
Point <- Point + tmp.point * tmp.weight
}
Points.unif[ii] <- Point
rs <- matrix(0,n,n.r/n)
for(i in 1:n){
rs[i,] <- rgamma(n.r/n,shape=sh[i],scale=sc[i])
}
Point <- 0
for(i in 1:length(rs)){
Point <- Point + dnorm(rs[i],M,S) * 1/length(rs)
}
Points.gamma[ii] <- Point
rs <- rnorm(n.r,mm,ss)
Point <- 0
for(i in 1:length(rs)){
tmp.prob <- 0
for(j in 1:n){
tmp.prob <- tmp.prob + dgamma(rs[i],shape=sh[j],scale=sc[j])
}
tmp.prob <- tmp.prob/n
tmp.point <- dnorm(rs[i],M,S) * 1/length(rs)
tmp.weight <- tmp.prob/dnorm(rs[i],mm,ss)
Point <- Point + tmp.point * tmp.weight
}
Points.norm[ii] <- Point
}
Points.three <- cbind(Points.unif,Points.gamma,Points.norm)
boxplot(Points.three)
> apply(Points.three,2,mean)
Points.unif Points.gamma Points.norm
6.670910e-06 6.483460e-06 6.389868e-06
> apply(Points.three,2,sd)
Points.unif Points.gamma Points.norm
2.819655e-06 9.335974e-06 7.721868e-06
-
- 今度は、正規乱数の方がばらつきがわずかだが小さくなったことが解る
- このようなことを勘案して、どんな乱数を使うのがよいかを検討してサンプリングするのがインポータンス・サンプリング
- ただし、生起確率密度がわかっていて、そのあたりの得点が高いときには、正規確率密度から乱数を発生させられるものなら発生させればよい。しかしながら、生起確率密度はわかるけれど、そこからの乱数発生関数がないような場合には、近くて乱数発生関数があるような乱数で代用し、その代り、重みづけについては、密度の比を取るのがよい
- また、生起確率密度からの乱数発生ができる場合であっても、大きな得点を得られるところの密度が低いと、乱数が無駄になる。そんなときには大きな得点が得られるところに高密度な別の乱数発生分布を設定して、それとの比を重みにするのがよい。インポータンスサンプリングの主の使い道はこれだろう(か？)
- それ以外の場合で、「得点を表した乱数」を発生することが良いような場合であっても、重みを付けるには、本当の生起確率密度が計算できなければ比が出せないことは同じ
- さて。Exponential tilting
- これは、生起確率分布が指数型分布族に属していて、乱数発生に適切な分布も指数型分布族に属しているときには、両者を指数型分布族表現してやり、その比(それがインポータンス・サンプリングの重み)が指数表現になっていて簡単だよ、という話
- 参考：こちらとか、こちらとか、こちらとか
