pquan<-function(pn,p,n){  
	#Returns Wilsons's Estimate  (Takes a sample proportion, 
		#P-value and sample size and returns a sample quantile)
	value<-NULL
	for(i in 1:length(p)){ 
		value[i]=(pn+qnorm(p[i])^2/(2*n)+
			(qnorm(p[i])/sqrt(n))*sqrt(pn*(1-pn)+qnorm(p[i])^2/(4*n)))/
				(1+(qnorm(p[i])^2/n))	
	}
	return(value)
}

pquanprox<-function(pn,p,n){
	value<-NULL
	for(i in 1:length(p)) value[i]=
		pn+(qnorm(p[i])/sqrt(n))*sqrt(pn*(1-pn))
	return(value)
}


betaP<-function(Qp,p,n,a=.5,b=.5){
	value<-NULL
	k<-Qp*n
	for(i in 1:length(p)) value=c(value,qbeta(p[i],k+a,n-k+b))
	return(value)
}

reliP<-function(Qp,p,a,b,n){
val<-n*(Qp*log(p/Qp)+(1-Qp)*log((1-p)/(1-Qp)))
return(exp(val))
}



op<-par(mfrow=c(1,1), mar=c(5,5,3,1),oma=c(0,0,0,0))

procbernoulli<-function(Qp,n,a=.5,b=.5,title="data"){
	P<-seq(1/(2*n),1-(1/(2*n)),length=1002)
	P<-P[2:1001]
	plot(P,pquan(Qp,P,n),main=title,,ylab="Quantile(p)",xlim=c(-.01,1),
		ylim=c(0,1.2),axes=F,col="steelblue4",type="l",lty=2)
		 #Wilson's estimate
	axis(1,seq(0,1,by=.1))
	axis(2,seq(0,1,by=.1))
	bp<-betaP(Qp,P,n,a,b)
	bn<-length(bp)
	px<-seq((1/(2*n)),(1-1/(2*n)),length=bn+2)
	px<-px[2:(bn+1)]
	points(px,bp,type="l",col="orangered3" ,lty=1)  #beta posterior 
	legend(0,1.2,c("Wilson's Estimate","beta-posterior"),
	lty=c(2,1),col=c("steelblue4","orangered3"),bty="n")
	text(.1,1.035,"beta-prior",pos=4)
	points(c(.085,.06,.035),c(1.035,1.035,1.035),pch=20,cex=.8)
	text(.1,1.2,"Relative Likelihood",pos=4)
	points(c(.085,.06,.035),c(1.2,1.2,1.2),pch=1,cex=.8)
	text(.65,1.15,"P =",pos=2,cex=1.25)
	text(.62,1.15,Qp,pos=4,cex=1.25)
	text(.85,1.15,"n =",pos=2,cex=1.25)
	text(.82,1.15,n,pos=4,cex=1.25)
	text(.65,1.05,"a =",pos=2,cex=1.25)
	text(.62,1.05,a,pos=4,cex=1.25)
	text(.85,1.05,"b =",pos=2,cex=1.25)
	text(.82,1.05,b,pos=4,cex=1.25)
	points(c(-.02,-.02),c(pquan(Qp,P[25],n),pquan(Qp,P[975],n)),
		col="steelblue4",type="o",pch=20,lty=2)
	points(c(-.02,.975),c(pquan(Qp,P[975],n),pquan(Qp,P[975],n)),
		col="steelblue4",type="o",pch=20,lty=2)
	#points(c(-.02,.025),c(pquan(Qp,P[25],n),pquan(Qp,P[25],n)),
		#col="steelblue4",type="o",pch=20,lty=2)
	points(c(-.01,-.01),c(betaP(Qp,P[25],n,a,b),betaP(Qp,P[975],n,a,b)),
		col="orangered3",type="o",pch=20)
	points(c(-.01,.975),c(betaP(Qp,P[975],n,a,b),betaP(Qp,P[975],n,a,b)),
		col="orangered3",type="o",pch=20)
	points(c(-.01,.025),c(betaP(Qp,P[25],n,a,b),betaP(Qp,P[25],n,a,b)),
		col="orangered3",type="o",pch=20)

	
	x<-qbeta(P,a,b)
	ind<-seq(1,1000,by=15)
	x=x[ind]
	P=P[ind]
	points(P,x,pch=19,cex=.6)#beta prior
	abline(h=Qp,lty=3,col="orchid3"  )

	P<-seq(1/(2*n),1-(1/(2*n)),length=302)
	P<-P[2:301]
	
	points(reliP(Qp,P,a,b,n),P,pch=1,cex=.8) #relative likelihood
}

procbernoulli(.2,n=36)