
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]))
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))
pp.sign <- sign(p1-p2)
diff.pp.sign <- diff(pp.sign)
diff.pp.sign2 <- diff.pp.sign[2:(length(diff.pp.sign)-1)]
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
}
out.optim <- optim(p.seq[tmp2],fx,lower=p.seq[tmp2],upper=p.seq[tmp3],method="L-BFGS-B")
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)