Perverting Burris-Hoff confidence intervals for fun and profit.

hoffci.R

## https://doi.org/10.1093/jssam/smz010

make.s<-function(mu, sigma,tau2,alpha){
	g<-make.g(alpha)
	invg<-make.ginv(alpha)
	
	function(theta) invg(2*sigma*(theta-mu)/tau2)
}

make.g<-function(alpha){
	function(omega)	qnorm(alpha*omega)-qnorm(alpha*(1-omega))
}
make.ginv<-function(alpha){
	g<-function(omega)	qnorm(alpha*omega)-qnorm(alpha*(1-omega))
	function(x){
		f<-function(omega) g(omega)-x
		if (x<0) {
	 		uniroot(f,c(pnorm(x/alpha),0.5))$root
		} else if (x==0) {
			0.5
		} else 
			uniroot(f,c(0.5,pnorm(x/alpha)))$root
	}
}

ci<-function(y, sigma, s,alpha){

	f_low<-function(theta) theta-(y+sigma*qnorm(alpha*(1-s(theta))))
	f_up<-function(theta) (y+sigma*qnorm(1-alpha*s(theta)))-theta

	low0<-y+qnorm(alpha/2)*sigma
	low1<-low0
	low<-if (f_low(low1)>0){
		while(f_low(low1)>0) low1<- y+(low1-y)*2
		uniroot(f_low,c(low1,low0))$root
	} else{
		uniroot(f_low,c(low0,y))$root
	}
	hi0<-y+qnorm(1-alpha/2)*sigma
	hi1<-hi0
	high<-if (f_up(hi1)>0){
		while(f_up(hi1)>0) hi1<- y+(hi1-y)*2
		uniroot(f_up,c(hi0,hi1))$root
	} else{
		uniroot(f_up,c(y,hi0))$root
	}
	
	c(low,high)
}

set.seed(as.Date("2019-6-7"))
s<-make.s(1,1,1,.1)
rr<-replicate(200,{a<-rnorm(1);c(a,ci(a,1, s, 0.1))})
cover<-(rr[2,]<0) & (rr[3,]>0)
plot(1:200, rr[1,],ylim=range(rr),pch=19,col=ifelse(cover,"grey","red"))
segments(1:200,y0=rr[2,],y1=rr[3,],col=ifelse(cover,"grey","red"))
abline(h=0,lty=2)