pquan<-function(pn,p,n){
	value<-NULL
	for(i in 1:length(p)){
	value=c(value,(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)
}

procpro<-function(s1,n1,s2,n2,title="data"){
	plot(1:150,seq(1,100,length=150),axes=F,main=title,type="n",
		xlab="",ylab="")
	points(c(50,50,150,150,50),c(0,100,100,0,0),type="l")
	points(c(100,100),c(0,100),type="l")
	points(c(50,150),c(50,50),type="l")
	text(75,75,s1,cex=5)
	text(75,25,s2,cex=5)
	text(125,75,n1-s1,cex=5)
	text(125,25,n2-s2,cex=5)
	text(25,75,"counts 1",cex=2)
	text(25,25,"counts 2",cex=2)
	
ck <- readline('strike "Enter" to continue')
	op<-par(mfrow=c(1,1),mar=c(5,4,4,1))
	ph1<-s1/n1
	ph2<-s2/n2
	P<-seq(0,1,length=102)
	P=P[2:101]
	p1<-NULL
	p2<-NULL
	for(i in 1:100){
		if(P[i]>(1/(2*n1))&&P[i]<1-(1/(2*n1))){
			p1=c(p1,pquan(ph1,P[i],n1))
		}
		if(P[i]>(1/(2*n2))&&P[i]<1-(1/(2*n2))){
			p2=c(p2,pquan(ph2,P[i],n2))
		}
	}
	lp1<-length(p1)
	lp2<-length(p2)
	P<-seq(1/(2*n1),1-(1/(2*n1)),length=lp1+2)
	P<-P[2:(lp1+1)]
	plot(P,p1,main=title,ylim=c(min( c(p1,p2) ),max(c(p1,p2))),
		pch=1,ylab="p")
	P<-seq(1/(2*n2),1-(1/(2*n2)),length=lp2+2)
	P<-P[2:(lp2+1)]
	points(P,p2,pch=16)
	ph1<-s1/n1
	ph2<-s2/n2
	P<-seq(0,1,length=10002)
	P=P[2:10001]
	d<-NULL
	z<-NULL
	for(i in 1:10000){
		z=runif(2)
		if(z[1]>(1/(2*n1)) && z[1]<(1-(1/(2*n1)))){
			if(z[2]>1/(2*n2)&&z[2]<(1-(1/(2*n2)))){
				d=c(d,pquan(ph2,z[1],n2)-pquan(ph1,z[2],n1))
			}
		}
	}
	d=sort(d)
	dn<-length(d)
	n<-min(n1,n2)
	P<-seq(1/(2*n),1-(1/(2*n)),length=dn+2)
	P<-P[2:(dn+1)]
	points(P,d,type="l")
	M<-max(p1,p2)
	points(c(.01,.02,.03,.04),rep(M,4))
	text(.05,M,expression(p[1]),pos=4)
	points(c(.01,.02,.03,.04),rep(M,4)-.025,pch=19)
	text(.05,M-.025,expression(p[2]),pos=4)
	points(c(.01,.04),rep(M,2)-.05,type="l")
	text(.05,M-.05,expression(p[2]-p[1]),pos=4)

ck <- readline('strike "Enter" to continue')
op<-par(mfrow=c(1,2),mar=c(1,1,4,1))
plot(1:100,1:100,type="n",axes=F,main="Quantile Points of p1",sub=title,
	xlab="",ylab="")
pm<-c(.5,.25,.1,.05,.025,.01,.005)
q1<-round(c(pquan(ph1,pm,n1)),4)
q2<-round(c(pquan(ph1,1-pm,n1)),4)
if(length(q1)<6)q1[6]="-----"
if(length(q1)<7)q1[7]="-----"
if(length(q2)<6)q2[6]="-----"
if(length(q2)<7)q2[7]="-----"
text(rep(15,5),seq(90,10,length=8),c("p",pm),pos=2)
text(rep(30,5),seq(90,10,length=8),c("Q(p)",q1),pos=4)
text(rep(80,5),seq(90,10,length=8),c("Q(1-p)",q2),pos=4)
 points(c(23,23),c(0,100),type="l")
abline(85,0)

plot(1:100,1:100,type="n",axes=F,main="Quantile Points of p2",sub=title,
	xlab="",ylab="")
pm<-c(.5,.25,.1,.05,.025,.01,.005)
q1<-round(c(pquan(ph2,pm,n2)),4)
q2<-round(c(pquan(ph2,1-pm,n2)),4)
if(length(q1)<6)q1[6]="-----"
if(length(q1)<7)q1[7]="-----"
if(length(q2)<6)q2[6]="-----"
if(length(q2)<7)q2[7]="-----"
text(rep(15,5),seq(90,10,length=8),c("p",pm),pos=2)
text(rep(30,5),seq(90,10,length=8),c("Q(p)",q1),pos=4)
text(rep(80,5),seq(90,10,length=8),c("Q(1-p)",q2),pos=4)
 points(c(23,23),c(0,100),type="l")
abline(85,0)

ck <- readline('strike "Enter" to continue')
op<-par(mfrow=c(1,1),mar=c(0,0,3,0))

plot(1:100,1:100,type="n",axes=F,main="Quantile Points p2-p1",sub=title,
	xlab="",ylab="")
pm<-c(.5,.25,.1,.05,.025,.01,.005)
q1<-c(d[round(.5*dn,0)],d[round(.25*dn,0)],d[round(.1*dn,0)],
	d[round(.05*dn,0)],d[round(.025*dn,0)],d[round(.01*dn,0)],
	d[round(.005*dn,0)])
q2<-c(d[round(.5*dn,0)],d[round(.75*dn,0)],d[round(.9*dn,0)],
	d[round(.95*dn,0)],d[round(.975*dn,0)],d[round(.99*dn,0)],
	d[round(.995*dn,0)])
text(rep(15,5),seq(90,10,length=8),c("p",pm),pos=2)
text(rep(30,5),seq(90,10,length=8),c("Q(p)",round(q1,4)),pos=4)
text(rep(80,5),seq(90,10,length=8),c("Q(1-p)",round(q2,4)),pos=4)
 points(c(23,23),c(0,100),type="l")
abline(85,0)

}

#op<-par(mfrow=c(1,1),mar=c(5,4,4,1))
#procpro(2,38,12,32,"stomach ulcers")

procodd<-function(s1,n1,s2,n2,title="data"){
	plot(1:150,seq(1,100,length=150),axes=F,main=title,type="n",xlab="",ylab="")
	points(c(50,50,150,150,50),c(0,100,100,0,0),type="l")
	points(c(100,100),c(0,100),type="l")
	points(c(50,150),c(50,50),type="l")
	text(75,75,s1,cex=5)
	text(75,25,s2,cex=5)
	text(125,75,n1-s1,cex=5)
	text(125,25,n2-s2,cex=5)
	text(25,75,"counts 1",cex=2)
	text(25,25,"counts 2",cex=2)
	
ck <- readline('strike "Enter" to continue')
	op<-par(mfrow=c(1,1),mar=c(5,4,4,1))
	ph1<-s1/n1
	ph2<-s2/n2
	n<-min(n1,n2)
	P<-seq(1/(2*n),(1-(1/(2*n))),length=102)
	P=P[2:101]
	y1<-NULL
	y2<-NULL
	y1s<-NULL
	y2s<-NULL
	for(i in 1:100){
		if(P[i]>(1/(2*n))&&P[i]<(1-(1/(2*n)))){
			y1s=pquan(ph1,P[i],n1)
			y1s=log(y1s/(1-y1s))
			y1=c(y1,y1s)
			y2s=pquan(ph2,P[i],n2)
			y2s=log(y2s/(1-y2s))
			y2=c(y2,y2s)
			}
		}
	plot(P,y1,main=title,ylab="theta",ylim=c( min(y1,y2) ,max(y1,y2)),
		pch=1)
	points(P,y2,pch=19)
	P<-seq(1/(2*n),(1-1/(2*n)),length=10002)
	P=P[2:10001]
	d<-NULL
	z<-NULL
	q1<-NULL
	q2<-NULL
	qd1<-NULL
	qd2<-NULL
	for(i in 1:10000){
		z=runif(2)
		if(z[1]>(1/(2*n1))&&z[1]<(1-(1/(2*n1)))){
			if(z[2]>(1/(2*n2))&&z[2]<(1-(1/(2*n2)))){
				q1=pquan(ph1,z[1],n1)
				q2=pquan(ph2,z[2],n2)
				d=c(d,log(q1/(1-q1))-log(q2/(1-q2)))
		}	}		
	}
	d<-sort(d)
	dn<-length(d)
	P<-seq(1/(2*n),(1-(1/(2*n))),length=dn+2)
	P=P[2:(dn+1)]
	points(P,d,type="l")
	points(c(.9,.91,.92,.93),c(-3.6,-3.6,-3.6,-3.6))
	text(.62,-3.6,"log_odds p1",pos=4)
	points(c(.9,.91,.92,.93),c(-3.9,-3.9,-3.9,-3.9),pch=19)
	text(.62,-3.9,"log_odds p2",pos=4)
	points(c(.895,.935),c(-4.2,-4.2),type="l")
	text(.62,-4.2,"log_odds p2-p1",pos=4)
	print(d[round(.025*dn,0)])
	print(d[round(.975*dn,0)])

ck <- readline('strike "Enter" to continue')
op<-par(mfrow=c(1,2),mar=c(1,1,4,1))
plot(1:100,1:100,type="n",axes=F,main="Quantile Points of p1",xlab="",
	ylab="")
pm<-c(.5,.25,.1,.05,.025,.01,.005)
q1<-pquan(ph1,pm,n1)
qd1<-pquan(ph1,1-pm,n1)
q2<-pquan(ph2,pm,n2)
qd2<-pquan(ph2,1-pm,n2)
for(i in 1:7){
	q1[i]=log(q1[i]/(1-q1[i]))
	q2[i]=log(q2[i]/(1-q2[i]))
	qd1[i]=log(qd1[i]/(1-qd1[i]))
	qd2[i]=log(qd2[i]/(1-qd2[i]))}
q1<-round(q1,4)
q2<-round(q2,4)
qd1<-round(qd1,4)
qd2<-round(qd1,4)
text(rep(15,5),seq(90,10,length=8),c("p",pm),pos=2)
text(rep(30,5),seq(90,10,length=8),c("Q(p)",q1),pos=4)
text(rep(80,5),seq(90,10,length=8),c("Q(1-p)",qd1),pos=4)
 points(c(23,23),c(0,100),type="l")
abline(85,0)

plot(1:100,1:100,type="n",axes=F,main=title,
	sub="quantile points(logodds of p1,p2)",xlab="",ylab="")
text(rep(15,5),seq(90,10,length=8),c("p",pm),pos=2)
text(rep(30,5),seq(90,10,length=8),c("Q(p)",q2),pos=4)
text(rep(80,5),seq(90,10,length=8),c("Q(1-p)",qd2),pos=4)
 points(c(23,23),c(0,100),type="l")
abline(85,0)

ck <- readline('strike "Enter" to continue')
op<-par(mfrow=c(1,1),mar=c(1,1,5,1))

plot(1:100,1:100,type="n",axes=F,main="Quantile Points log odds p1-p2",
	xlab="",ylab="")
pm<-c(.5,.25,.1,.05,.025,.01,.005)
q1<-c(d[round(.5*dn,0)],d[round(.25*dn,0)],d[round(.1*dn,0)],
	d[round(.05*dn,0)],d[round(.025*dn,0)],d[round(.01*dn,0)],
	d[round(.005*dn,0)])
q2<-c(d[round(.5*dn,0)],d[round(.75*dn,0)],d[round(.9*dn,0)],
	d[round(.95*dn,0)],d[round(.975*dn,0)],d[round(.99*dn,0)],
	d[round(.995*dn,0)])
text(rep(15,5),seq(90,10,length=8),c("p",pm),pos=2)
text(rep(30,5),seq(90,10,length=8),c("Q(p)",round(q1,4)),pos=4)
text(rep(80,5),seq(90,10,length=8),c("Q(1-p)",round(q2,4)),pos=4)
 points(c(23,23),c(0,100),type="l")
abline(85,0)


}


op<-par(mfrow=c(1,1),mar=c(5,4,4,1))
procpro(2,38,12,32,"stomach ulcers")

op<-par(mfrow=c(1,1),mar=c(5,4,4,1))
procodd(2,38,12,32,"quantile plot of log odds p1,p2,and p2-p1")