op<-par(mfrow=c(1,1),mar=.1+c(6,4,4,1))

procnormal<-function(ybar,mu0=1,sig0=100,sig=1,n=100,title="data"){
	P<-1:100
	mug<-seq(1,3,length=300)
	P<-(P-.5)/100
	q0<-mu0+sig0*qnorm(P)
	plot(P,q0,xlab="P-value",ylab="means",axes=F,main=title,
		sub="Confidence Quantiles",ylim=c(1,3.5),
			xlim=c(-.02,1),pch=20)#prior
	axis(1,seq(0,1,by=.1))
	axis(2,seq(1,3,by=.5))
	legend(.05,3.5,c("Confidence Quantile","posterior"),
		lty=c(2,1),bty="n")
	text(.6,3.4,"prior",pos=4)
	points(c(.55,.57,.59),c(3.4,3.4,3.4),pch=20)
	text(.6,3.3,"relative likelihood",pos=4)
	points(c(.55,.57,.59),c(3.3,3.3,3.3),pch=1)
	I0<-1/sig0^2
	In<-n/sig^2
	I1<-I0+In
	sig1<-sqrt(1/I1)
	mu1<-(I0/I1)*mu0+ybar*(In/I1)
	q1<-mu1+sig1*qnorm(P)
	points(P,q1,type="l")#posterior
	abline(ybar,0,lty=2,col="springgreen4")
	p<-exp(-.5*(n/sig^2)*(mug-ybar)^2)

	points(p,mug,pch=1)#relative likelihood
	y2<-mu1+sig1*qnorm(.025)
	y2[2]=mu1+sig1*qnorm(.975)
	points(c(-.01,-.01),y2,type="o",pch=9)
	y<-ybar+sqrt(2*(-log(.135))*sig^2/n)
	y[2]=ybar-sqrt(2*(-log(.135))*sig^2/n)
	points(P,qnorm(P,ybar,sig/sqrt(n)),type="l",lty=2)
			#Confidence Quantile 
	points(c(-.03,-.03),y,type="o",pch=20,col="orangered4")
	points(c(.5,.5),c(0,3),type="l",lty=3,col="khaki4")
	points(c(-.03,.135),c(y[1],y[1]),type="l",col="red")
	points(c(-.03,.135),c(y[2],y[2]),type="l",col="red")
	points(c(-.01,.025),c(y2[1],y2[1]),type="l")
	points(c(-.01,.975),c(y2[2],y2[2]),type="l")
	points(c(-.05,.025),c(qnorm(.025,ybar,sig/sqrt(n)),
		qnorm(.025,ybar,sig/sqrt(n))),type="l",lty=2)
	points(c(-.05,.975),c(qnorm(.975,ybar,sig/sqrt(n)),
		qnorm(.975,ybar,sig/sqrt(n))),type="l",lty=2)
	points(c(-.05,-.05),c(qnorm(.975,ybar,sig/sqrt(n)),
		qnorm(.025,ybar,sig/sqrt(n))),type="o",lty=2)
	text(.1,3,expression(mu[0]),cex=1.25)
	text(.14,3,c(" = "),cex=1.25)
	text(.14,3.01,mu0,pos=4)
	text(.35,3,expression(sigma[0]),cex=1.25)
	text(.39,3,c(" = "),cex=1.25)
	text(.39,3.01,sig0,pos=4)
	text(.6,3,expression(n),cex=1.25)
	text(.64,3,c(" = "),cex=1.25)
	text(.64,3.01,n,pos=4)
	text(.1,2.85,expression(mu[1]),cex=1.25)
	text(.14,2.85,c(" = "),cex=1.25)
	text(.14,2.851,rd(mu1),pos=4)
	text(.35,2.85,expression(sigma[1]),cex=1.25)
	text(.39,2.85,c(" = "),cex=1.25)
	text(.39,2.851,rd(sig1),pos=4)
	text(.60,2.85,expression(bar(y)),cex=1.25)
	text(.64,2.85,c(" = "),cex=1.25)
	text(.64,2.851,ybar,pos=4)
	text(.8,3,expression(sigma),cex=1.25)
	text(.84,3,c(" = "),cex=1.25)
	text(.84,3,sig,pos=4)
}
procnormal(2,1.5,.2,1,n=100,"data")