２つのベータ分布

• ベータ関数とベータ分布についての基本
  • まずはWikiへのリンク
    • ベータ関数
    • 不完全ベータ関数
    • ベータ関数英語版(英語版では同じページに
  • ベータ分布の確率密度関数と累積密度関数は
    • 確率密度関数
    • 累積密度関数
      • ただしは不完全ベータ関数
      • はベータ関数で
      • は正則ベータ関数で
      • が整数のとき、

a <-3
b <- 4
beta(a,b)
pbeta(0.3,a,b) # 正則ベータ関数(ベータ分布の累積密度関数)
pbeta(0.3,a,b) * beta(a,b) # 不完全ベータ関数(B(x=0.3;a=3,b=4)

• ベータ分布
  • 真偽の事象で真偽がそれぞれ回観測されたときの、真の生起確率の分布としてベータ分布を共役事前分布としてとることがある

plot(dbeta(seq(from=0,to=1,length=100),a,b))

• ２つの相互に独立なベータ分布の積(複合分布)

X <- matrix(c(3,4,7,9),2,2)
x <- seq(from=0,to=1,length=100)
p1 <- dbeta(x,X[1,1]+1,X[1,2]+1)
p2 <- dbeta(x,X[2,1]+1,X[2,2]+1)
p12 <- outer(p1,p2,"*")
image(p12,xlim=c(0,1),ylim=c(0,1))
par(new=TRUE)
contour(p12,add=TRUE)

• ２ベータ分布の勝敗分布
  • ここで２分布の勝負を考える
    • 片方の分布からの値が、もう片方の分布からの値より大きくなる確率は、この複合分布において、対角線の片側の部分の積分

abline(0,1)

• • 式で表すと
  • 
  • 
  • 
  • 上に記した、正則ベータ関数とかを使って、式変形していく
  • 
  • 
  • 
  • なので
  • 
  • 
• Rで確かめてみよう
  • まずは、x,y,z,w４数を引数としてPを計算する関数を２つ
    • 式にべったりのそれと、の影響の有無でまとめ直して、対数が取れるところは対数をとった２つを
  • 注。引数がx=(x1,x2,x3,x4)の場合とx=(x4,x3,x2,x1)の場合とは、複合確率密度関数の積分の仕方の軸を変更してはいるが、同じ積分なので値は同じ。x1,x4の大小が、の項の数を決めているので、少ない方を採用するのが軽そうに思えるのでそうしておく

Decision_beta <- function(x){
	if(x[1] < x[4]){
		x <- x[4:1]
	}
	ret <- 0
	for(j in x[3]:(x[3]+x[4]-1)){
		ret <- ret + beta(x[1]+j,x[2]+x[3]+x[4]-1-j)/beta(j+1,x[3]+x[4]-j)/(x[3]+x[4])
	}
	return(ret/beta(x[1],x[2]))
}
Decision_beta.2 <- function(x){
	if(x[1] < x[4]){
		x <- x[4:1]
	}
	ret <- 0
	common <- -log(x[3]+x[4])+lgamma(x[1]+x[2])+lgamma(x[3]+x[4]+1)-lgamma(x[1])-lgamma(x[2])-lgamma(sum(x)-1)
	for(j in x[3]:(x[3]+x[4]-1)){
		tmp <- lgamma(x[1]+j)+lgamma(sum(x[2:4])-1-j)-lgamma(j+1)-lgamma(x[3]+x[4]-j)
		ret <- ret + exp(tmp+common)
	}
	return(ret)
}

• • ２つのベータ分布が真偽観測から作られることを考慮して、観測データを与えるバージョン。上の作りに合わせて異なる関数(結果は同じ)を３つ。Decision_Dice_beta.2.2()が一番よいだろう、多分。

Decision_Dice_beta.2 <- function(X){
	x <- X[1,1]
	y <- X[1,2]
	z <- X[2,1]
	w <- X[2,2]
	Decision_beta(c(x,y,z,w)+1)
}

Decision_Dice_beta.2.2 <- function(X){
	x <- X[1,1]
	y <- X[1,2]
	z <- X[2,1]
	w <- X[2,2]
	Decision_beta.2(c(x,y,z,w)+1)
}

Decision_Dice_beta <- function(X){
	x <- X[1,1]
	y <- X[1,2]
	z <- X[2,1]
	w <- X[2,2]
	ret <- 0
	for(j in (z+1):(z+w+1)){
		ret <- ret + beta(j+x+1,z+w+y-j+2)/beta(j+1,z+w+2-j)/(z+w+2)
	}
	return(ret/beta(x+1,y+1))
}

• Stirlingの近似は、将来的に大きな値に対応する必要が出たときに入用なので、Wikiのベータ関数から引用しておく
  • 

２つのベータ分布の勝敗確率関数のできを確かめる

• 前の記事で式変形したり関数を作ったりしたけれど、さて、これがあっているもんだか不安
• 複数の関数の返り値が一致しても、それは、式変形した結果得られたものをさらに変形しているだけなので、何の役にも立たない
• ベータ分布からの乱数を発生させて、勝負をさせて、そのうえで勝敗確率がいい感じになるかどうかを比べて確かめることにする
• ベータ分布からの乱数発生法の方は、２項から多項への展開を見越して、ディリクレ乱数としておく

Decision_Dice <- function(X,v=c(0,1),n=10000){
	R <- list()
	V <- matrix(0,n,length(X[,1]))
	for(i in 1:length(X[,1])){
		R[[i]] <- rdirichlet(n,X[i,]+1)
		V[,i] <- v %*% t(R[[i]])
	}
	S <- tabulate(apply(V,1,which.max),length(X[,1]))
	S/sum(S)
}

• モンテカルロ比較で、大丈夫そう

n.iter <- 1000

out1<-out2<-rep(0,n.iter)
for(i in 1:n.iter){
	X <- matrix(sample(0:100,4),2,2)
	#print("---")
	out1[i] <- Decision_Dice(X,c(1,0),n=10000)[1]
	out2[i] <- Decision_Dice_beta.2.2(X)
}
plot(out1,out2)

• さて、多項間の決断
  • 上のようにベータ分布のn次複合分布をこちらにあるように、n次元立方体のn等分錐に関して積分したいわけだが、正則ベータ関数の展開がもかもかと大きくなってしまって(もしかしたらうまく行くのかもしれないけれど、ぱっと見た感じ、うまく片付けられそうもないので)、モンテカルロで代用するのが賢明なのではないかと思われる