ベータ分布の交叉点

• ベータ分布が二つあるとする
• その累積分布を考える。ベータ分布の累積分布は、regularized incomplete beta 関数
• ２つのベータ分布の累積分布関数の交点を計算してみる

# x1は一つ目のベータ分布を決める、二つの成否回数ベクトル
# x2は二つ目のベータ分布を決める、二つの成否回数ベクトル
# p.seqはベータ分布の確率・累積確率を計算するサポートの値
da.description <- function(x1,x2,p.seq=NULL){
	if(is.null(p.seq)){
		p.seq <- seq(from=0,to=1,length=101)
	}
# 確率密度分布と累積分布との値を計算	
	d1 <- dbeta(p.seq,x1[1]+1,x1[2]+1)
	d2 <- dbeta(p.seq,x2[1]+1,x2[2]+1)
	p1 <- pbeta(p.seq,x1[1]+1,x1[2]+1)
	p2 <- pbeta(p.seq,x2[1]+1,x2[2]+1)
# サンプルの成功割合
	sample.success <- c(x1[1]/(x1[1]+x1[2]),x2[1]/(x2[1]+x2[2]))
# 観測がないときは、成功割合を0.5とみなすことにしている
	if(x1[1]+x1[2]==0){
		sample.success[1] <- 0.5
	}
	if(x2[1]+x2[2]==0){
		sample.success[2] <- 0.5
	}
# ベータ分布の期待値
	exp.success <- c((x1[1]+1)/(sum(x1)+2),(x2[1]+1)/(sum(x2)+2))
# ２つの累積分布の真の交点付近のp.seq値を交点の前後の値として取り出したい
# まず、２つの分布の大小を確認
	pp.sign <- sign(p1-p2)
# 大小関係が変化する点を確認
	diff.pp.sign <- diff(pp.sign)
	diff.pp.sign2 <- diff.pp.sign[2:(length(diff.pp.sign)-1)]
# そのp.seqの番地を取り出す
	tmp <- which(diff.pp.sign2!=0)[1]

	tmp2 <- tmp+1
	tmp3 <- tmp+2
# 前後値を探索範囲として、交点座標の最適化探索をする
# そのために最小化したい関数を定義
	fx <- function(x){
		(pbeta(x,x1[1]+1,x1[2]+1)-pbeta(x,x2[1]+1,x2[2]+1))^2
	}
# optim()に処理させる
	out.optim <- optim(p.seq[tmp2],fx,lower=p.seq[tmp2],upper=p.seq[tmp3],method="L-BFGS-B")
# optim()関数の出力とともに、色々を返す
	return(list(p.seq=p.seq,d=cbind(d1,d2),p=cbind(p1,p2),sample.success=sample.success,exp.success=exp.success,p.cross.optim=out.optim))
}

• 使ってみる

x1 <- c(5,3)
x2 <- c(25,25)
out.1 <- da.description(x1,x2)
out.1$p.cross.optim$par

matplot(out.1$p.seq,out.1$d,type="l")
abline(v=out.1$p.cross.optim$par,col=3)

matplot(out.1$p.seq,out.1$p,type="l")
abline(v=out.1$p.cross.optim$par,col=3)