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

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

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等分錐に関して積分したいわけだが、正則ベータ関数の展開がもかもかと大きくなってしまって(もしかしたらうまく行くのかもしれないけれど、ぱっと見た感じ、うまく片付けられそうもないので)、モンテカルロで代用するのが賢明なのではないかと思われる