# March 9, 2017 min 1/a - 1/b # March 8, 2017 # Compare widths of confidence intervals chisq2 = function(df=5,alpha=.05) { add.plt = function(df,a,b) { xx = seq(a,b,,501); yy = dchisq(xx,df) xx=c(a,xx,b); yy=c(0,yy,0); polygon(xx,yy,col=2) } crit = function(x,df,alpha) { a=exp(x[1]); b=exp(x[2]) fa=dchisq(a,df); fb=dchisq(b,df) pa=pchisq(a,df); pb=pchisq(b,df) val = (1/a-1/b)^2 + 1000*(a^2*fa-b^2*fb)^2 + 1000*(pb-pa-1+alpha)^2 if(br) {browser()}; return(val) } par(mfrow=c(3,1),mar=c(0,3,0,0),oma=c(5,0,5,0)) xx = seq(0,qchisq(.995,df),,501); fx=dchisq(xx,df) # equal tail a = qchisq(alpha/2,df); b = qchisq(1-alpha/2,df) plot(xx,fx,type="l",xaxt="n"); add.plt(df,a,b) abline(h=0); text(xx[350],max(fx)/1.5,paste("1/a-1/b =", round(1/a-1/b,3)),cex=1.5) axis(1,outer=T); title(paste("chisq df = ", df, " alpha =",alpha),outer=T,cex.main=2) # right tail a1 = 0; b1 = qchisq(1-alpha,df) plot(xx,fx,type="l",xaxt="n",yaxt="n"); add.plt(df,a1,b1) abline(h=0); text(xx[350],max(fx)/1.5,paste("1/a-1/b =", round(1/a1-1/b1,3)),cex=1.5) # equal height ans=nlminb(c(log(a),log(b)),crit,df=df,alpha=alpha) a2 = exp(ans$par[1]); b2 = exp(ans$par[2]) browser() plot(xx,fx,type="l",xaxt="n"); add.plt(df,a2,b2) abline(h=0); text(xx[350],max(fx)/1.5,paste("1/a-1/b =", round(1/a2-1/b2,3)),cex=1.5) abline(h=dchisq(a2,df),lty=3) browser() } chisq2( 5,0.10) chisq2(10,0.10) chisq2(20,0.10) chisq2(40,0.10) chisq2( 3,0.10)