# routine to compute 4 prob intervals and 4 conf intervals for the
# non-central chi distribution.  Taken from  
# Kent, J.T. and Hainsworth, T.J. (1995) Confidence intervals for the
# noncentral chi-squared distribution. J Stat Planning Inf 46, 147-159.

# Consider y^2 ~ \chi^2_p(\lambda^2)

# Given \lambda, we can find a (1-\alpha) prob interval for y, or
# given y, we can find a (1-\alpha) confidence interval for \lambda.

# Four methods for each are carried out by the function below.
# The results are summarized in an 8 x 3 matrix, where the first column
# gives the lower endpoint, the second column gives the upper endpoint,
# and the third column is a label for the type of interval.


  F=function(u,p,lambda) pchisq(u^2,p,lambda^2) # cdf of ncc
  Finv=function(prob,p,lambda) sqrt(qchisq(prob,p,lambda^2)) # quantile
  lambound=function(y,p,alpha) {
    lambda=1; obj=1
    while(obj>0) {
      lambda=2*lambda
      obj=F(y,p,lambda)-alpha
    }
    lambda
  }
lamfind=function(y,p,alpha) {
  if(alpha<1e-4 | F(y,p,0)<alpha) break("bad lamfind")
  lbig=lambound(y,p,alpha)
  f=function(lambda) F(y,p,lambda) - alpha
  lambda=uniroot(f,c(0,lbig))$root
}

  gb=function(u,nu,lambda) {
    # noncentral chi pdf wrt  Bessel base
    case1=(u>0 & lambda>0)
    case2=(u>0 & lambda==0)
    case3=(u==0)
    if(case1) pdf=lambda^(-nu)* exp(-.5*lambda^2)*besselI(lambda*u,nu)/
      sqrt(besselI(u^2,nu))
    if(case2) pdf=(u/2)^nu/(gamma(nu+1)*sqrt(besselI(u^2,nu)))
    if(case3) pdf=exp(-.5*lambda^2)/sqrt(gamma(nu+1))/(2^(nu/2))
    pdf
  }
  gr=function(u,nu,lambda) {
    # noncentral chi pdf wrt radial base
    ul=lambda*u
    if(ul>0) pdf=ul^(-nu)*besselI(ul,nu)*exp(-.5*(u^2+lambda^2))
    else pdf= (.5^nu)*exp(-.5*(u^2+lambda^2))/gamma(nu+1)
    pdf
  }
  # central prob int
ncc.pi.central=function(lambda,p,alpha=0.05) {
  int1=sqrt(qchisq(c(alpha/2,1-alpha/2),p,lambda^2))
  int1
}
  # MD Bessel prob int 
ncc.pi.mdb=function(lambda,p,alpha=0.05) {
  nu=(p-2)/2; u0l=Finv(1-alpha, p,lambda)
  if(gb(0,nu,lambda)>=gb(u0l,nu,lambda)) int2=c(0,u0l)
  else {
    f=function(cd) (log(gb(cd[1],nu,lambda)/gb(cd[2],nu,lambda)))^2 +
         (F(cd[2],p,lambda)-F(cd[1],p,lambda) - (1-alpha))^2
    int2=constrOptim(ncc.pi.central(lambda,p,alpha),f,grad=NULL,
      ui=diag(2),ci=c(0,0))$par
  }
  int2
}
  # MD radial prob int # special needed
ncc.pi.mdr=function(lambda,p,alpha=0.05) {
  nu=(p-2)/2; u0l=Finv(1-alpha,p, lambda)
  if(gr(0,nu,lambda)>=gr(u0l,nu,lambda)) int3=c(0,u0l)
  else {
    f=function(cd) (log(gr(cd[1],nu,lambda)/gr(cd[2],nu,lambda)))^2 +
         (F(cd[2],p,lambda)-F(cd[1],p,lambda) - (1-alpha))^2
    int3=constrOptim(ncc.pi.central(lambda,p,alpha),f,grad=NULL,
      ui=diag(2),ci=c(0,0))$par
  }
  int3
}
  # symmetric range prob int
ncc.pi.sr=function(lambda,p,alpha=0.05) {
  if(F(2*lambda,p,lambda)<=1-alpha) int4=c(0,Finv(1-alpha,p, lambda))
  else {
    f=function(b) F(lambda+b,p,lambda)-F(lambda-b,p,lambda)-(1-alpha)
    bb=uniroot(f,c(0,lambda))$root
    int4=c(lambda-bb,lambda+bb)
  }
  int4
}
  # central conf int
ncc.ci.central=function(y,p,alpha=0.05) {
  u1=Finv(alpha/2,p,0); u2=Finv(1-alpha/2,p,0)
  if(y<=u1) {ll=NaN; lu=NaN}
  else {
    lu=lamfind(y,p,alpha/2)
#    f=function(lambda) F(y,p,lambda)-alpha/2
#    lbig=lambound(y,p,alpha/2)
#    lu=uniroot(f,c(0,lbig))$root
    if(y <=u2) ll=0
    else ll=lamfind(y,p,1-alpha/2)
#    f=function(lambda) F(y,p,lambda)-(1-alpha/2)
#    lbig=lambound(y,p,1-alpha/2)
#    ll=uniroot(f,c(0,lbig))$root
  }
  int5=c(ll,lu)
  int5
}
  # MD Bessel conf int
ncc.ci.mdb=function(y,p,alpha=0.05) {
  u0=Finv(1-alpha,p,0); nu=(p-2)/2; zz=-qnorm(alpha/2)
  if(y<=u0) ll=0
  else {
    ly=lamfind(y,p,1-alpha)
#    f=function(lambda) F(y,p,lambda)-(1-alpha)
#    lbig=lambound(y,p,1-alpha)
#    ly=uniroot(f,c(0,lbig))$root
    g0=gb(0,nu,ly); gy=gb(y,nu,ly)
    if(g0>=gy) ll=ly
    else {
      f=function(cd) (log(gb(y,nu,cd[1])/gb(cd[2],nu,cd[1])))^2 +
           (F(y,p,cd[1])-F(cd[2],p,cd[1]) - (1-alpha))^2 # cd1=lam; cd2=cc
      start=c(y-zz,y-2*zz)
      if(start[1]<=0) start[1]=.1*(y+.1)
      if(start[2]<=0) start[2]=.05*(y+.1)
      ll=constrOptim(start,f,grad=NULL,ui=diag(2),ci=c(0,0))$par[1]
    }
  }
  f=function(cd) (log(gb(y,nu,cd[1])/gb(cd[2],nu,cd[1])))^2 +
       (F(cd[2],p,cd[1])-F(y,p,cd[1]) - (1-alpha))^2 # cd1=lam; cd2=d
  start=c(y+zz,y+2*zz)
  lu=constrOptim(start,f,grad=NULL,ui=diag(2),ci=c(0,0))$par[1]
  int6=c(ll,lu)
  int6
}
  # MD radial conf int
ncc.ci.mdr=function(y,p,alpha=0.05) {
  u0=Finv(1-alpha,p,0); nu=(p-2)/2; zz=-qnorm(alpha/2)
  if(y<=u0) ll=0
  else {
    ly=lamfind(y,p,1-alpha)
#    f=function(lambda) F(y,p,lambda)-(1-alpha)
#    lbig=lambound(y,p,1-alpha)
#    ly=uniroot(f,c(0,lbig))$root
    g0=gr(0,nu,ly); gy=gr(y,nu,ly)
    if(g0>=gy) ll=ly
    else {  
      f=function(cd) (log(gr(y,nu,cd[1])/gr(cd[2],nu,cd[1])))^2 +
           (F(y,p,cd[1])-F(cd[2],p,cd[1]) - (1-alpha))^2 # cd1=lam; cd2=cc
      start=c(y-zz,y-2*zz)
      if(start[1]<=0) start[1]=.1*(y+.1)
      if(start[2]<=0) start[2]=.05*(y+.1)      
      ll=constrOptim(start,f,grad=NULL,ui=diag(2),ci=c(0,0))$par[1]
    }
  }
  f=function(cd) (log(gr(y,nu,cd[1])/gr(cd[2],nu,cd[1])))^2 +
       (F(cd[2],p,cd[1])-F(y,p,cd[1]) - (1-alpha))^2 # cd1=lam; cd2=d
  start=c(y+zz,y+2*zz)
  lu=constrOptim(start,f,grad=NULL,ui=diag(2),ci=c(0,0))$par[1]
  int7=c(ll,lu)
  int7
}
  # symmetric range CI
ncc.ci.sr=function(y,p,alpha=0.05) {
  u0=Finv(1-alpha,p,0)  
  if(y<=u0) ll=0
  else {
    f=function(b) F(y,p,y-b)-F(max(y-2*b,0),p,y-b)-(1-alpha)
    bb=uniroot(f,c(0,y))$root
    ll=y-bb
  }
  f=function(b) F(y+2*b,p,y+b)-F(y,p,y+b)-(1-alpha)
  big=1; obj=-1
  while(obj<0) {big=2*big; obj=f(big)}
  bb=uniroot(f,c(0,big))$root
  lu=y+bb
  int8=c(ll,lu)
  int8
}
ncc.ints=function(yl,p,alpha=0.05) {
  int1=ncc.pi.central(yl,p,alpha)
  int2=ncc.pi.mdb(yl,p,alpha)
  int3=ncc.pi.mdr(yl,p,alpha)
  int4=ncc.pi.sr(yl,p,alpha)
  int5=ncc.ci.central(yl,p,alpha)
  int6=ncc.ci.mdb(yl,p,alpha)
  int7=ncc.ci.mdr(yl,p,alpha)
  int8=ncc.ci.sr(yl,p,alpha)
  ans=round(rbind(int1,int2,int3,int4,int5,int6,int7,int8),2)
  rownames(ans)=c("Central prob int","MD Bessel prob int",
                       "MD radial prob int","Symmetric range prob int",
                       "Central conf int","MD Bessel conf int",
                       "MD radial conf int","Symmetric range conf int")
  ans
}

test=function(lambda,p,alpha=0.05) {
  int1=ncc.pi.central(lambda,p,alpha)
  int5l=ncc.ci.central(int1[1],p,alpha); int5u=ncc.ci.central(int1[2],p,alpha)
  print(c(int1,int5l,int5u))
  int2=ncc.pi.mdb(lambda,p,alpha)
  int6l=ncc.ci.mdb(int2[1],p,alpha); int6u=ncc.ci.mdb(int2[2],p,alpha)
  print(c(int2,int6l,int6u)) 
  int3=ncc.pi.mdr(lambda,p,alpha)
  int7l=ncc.ci.mdr(int3[1],p,alpha); int7u=ncc.ci.mdr(int3[2],p,alpha)
  print(c(int3,int7l,int7u))
  int4=ncc.pi.sr(lambda,p,alpha)
  int8l=ncc.ci.sr(int4[1],p,alpha); int8u=ncc.ci.sr(int4[2],p,alpha)
  print(c(int4,int8l,int8u))
  NULL
}