18/11/10 09:12:19.51 BVAJndFJ.net
seq2pCI <- function(N,K,alpha=0.05){
vp=Vectorize(function(p)seqNp(N,K,p)-seqNp(N,K+1,p))
curve(vp(x),lwd=2,bty='l',xlab='Pr[head]',ylab=paste('Pr[max',K,'-head repetition]'))
abline(h=alpha,lty=3)
peak=optimize(vp,c(0,1),maximum=TRUE)$maximum
mean=integrate(function(x)x*vp(x),0,1)$value/integrate(function(x)vp(x),0,1)$value
lwr=uniroot(function(x,u0=alpha) vp(x)-u0,c(0.01,peak))$root
upr=uniroot(function(x,u0=alpha) vp(x)-u0,c(peak,0.99))$root
c(lower=lwr,mean=mean,mode=peak,upper=upr)
}
seq2pCI(100,5)