################################################################################
## max asymptotic MSE (psi)                                                     
#[Ru1] Higher Order Expansion for the MSE of M-estimators on shrinking neighborhoods
#[Ru2] Consequences of Higher Order Asymptotics for the MSE of M-estimators on Neighborhoods
# (c) P. Ruckdeschel May-2010                                                   
################################################################################


asMSE0<-function(r,b,v0)  ## [Ru1](3.9)
  {r^2*b^2+v0^2}

asMSE1<-function(n,r,b,v0,l2,v1,pm) ## [Ru1](3.10)
  {A1<-v0^2*(pm*(4*v1+3*l2)*b+1)+b^2+(2*b^2+pm*l2*b^3)*r^2
   asMSE0(r,b,v0)+A1*r/sqrt(n)}

asMSE2<-function(n,r,b,v0,l2,v1,pm,l3,v2,r1,r0) ## [Ru1](3.11)
  {A2<-v0^3*((l2+2*v1)*r0+2/3*r1)+v0^4*(3*v2+15/4*l2^2+l3+9*v1^2+12*v1*l2)
   A2<-A2+(v0^2*((3*v2+3*v1^2+15/2*l2^2+2*l3+12*v1*l2)*b^2+1+pm*(8*v1+6*l2)*b))*r^2
   A2<-A2+(pm*3*l2*b^3+5*b^2)*r^2+((5/4*l2^2+l3/3)*b^4+pm*3*l2*b^3+3*b^2)*r^4
   asMSE1(n,r,b,v0,l2,v1,pm)+A2/n}

##for Hampel-type ICs in Gaussian location model:


Vnorm<-function(c){Anorm(c)^2*(2*pnorm(c)-1-2*dnorm(c)*c+2*c^2*pnorm(-c))} #v0 [Ru1](Prop.3.4)
Anorm<-function(c){(2*pnorm(c)-1)^(-1)}        #A [Ru1](Prop.3.4)
l3norm<-function(c){2*Anorm(c)*c*dnorm(c)}     #l3 [Ru1](Prop.3.4)
v2norm<-function(c){p<-pnorm(c); d<-dnorm(c);  #v2 [Ru1](Prop.3.4)
                    Z<-6*p-4*p^2-2-2*c*d; 
                    N<-2*c^2*(1-p)+2*p-1-2*c*d; Z/N}
r1norm<-function(c){p<-pnorm(c); d<-dnorm(c); A<-Anorm(c); v0<-sqrt(Vnorm(c)); #rho1(Prop.3.4)
                    3*A^3*(1-2*p+2*c*d)/v0^3+3/v0}

asMSE0.norm<-function(r,c){  A<-Anorm(c); b=c*A; v0<-sqrt(Vnorm(c)); asMSE0(r,b,v0)}
asMSE1.norm<-function(n,r,c){A<-Anorm(c); b=c*A; v0<-sqrt(Vnorm(c)); asMSE1(n,r,b,v0,0,0,1)}
asMSE2.norm<-function(n,r,c){A<-Anorm(c); b=c*A; v0<-sqrt(Vnorm(c)); 
                             l3<-l3norm(c);v2<-v2norm(c);r1<-r1norm(c);
                             asMSE2(n,r,b,v0,0,0,1,l3,v2,r1,0)}

###optimal clipping heights: 

c0opt<-function(r)#  [Ru2](2.10)
#calculates c0
{fs<-function(k,r){-k+k*pnorm(k)+dnorm(k)-r^2/2*k}
 uniroot(fs,low=0,up=10,tol=10^-7,r=r)$root
}

c1opt<-function(n,r)#  [Ru2](4.5)
#calculates c1
{fs<-function(k,r,n){-k+k*pnorm(k)+dnorm(k)-r^2/2*k*(1+(r^2+1)/(r^2+r*sqrt(n)))}
 uniroot(fs,low=0,up=10,tol=10^-7,r=r,n=n)$root
}

c2opt<-function(n,r)# brute force
{optimize(asMSE2.norm,low=0,up=10,tol=10^-7,r=r,n=n)$minimum
}


### minimax radius:

rho<-function(r1,r0,ord,nn){  ## of all orders
#
if(ord==0) {c1<-c0opt(r1);c0<-c0opt(r0); 
            MSE<-function(nn,r,c)asMSE0.norm(r,c) }
if(ord==1) {c1<-c1opt(nn,r1);c0<-c1opt(nn,r0);
            MSE<-function(nn,r,c)asMSE1.norm(nn,r,c) }
if(ord==2) {c1<-c2opt(nn,r1);c0<-c2opt(nn,r0);
            MSE<-function(nn,r,c)asMSE1.norm(nn,r,c) }
MSE(nn,r0,c1)/MSE(nn,r0,c0)
}



getrmax<-function(r1,ord,nn){
     rma<-ifelse(ord==0,upp,sqrt(nn))
     max(rho(r1=r1,r=10^-6,ord,nn),rho(r1=r1,r=rma,ord,nn))
     }
getr0<-function(ord,nn){
   rma<-ifelse(ord==0,upp,sqrt(nn))
   erg<-optimize(getrmax,lower=10^-10,upper=rma,ord=ord,nn=nn)
   rho<-erg$objective; r0<-erg$minimum;
   if(ord==0) cr0<-c0opt(r0)
   if(ord==1) cr0<-c1opt(nn,r0)
   if(ord==2) cr0<-c2opt(nn,r0)
   c(rho=rho,r0=r0,cr0=cr0)
}


rhogam<-function(r1,r00,gam0,ord,nn){
      rma<-ifelse(ord==0,r00*gam0,min(sqrt(nn),r00*gam0))
      max(rho(r1=r1,r0=r00/gam0,ord,nn),rho(r1=r1,r0=rma,ord,nn))
      }

rhogamo<-function(r,gam0,ord=ord,nn=nn){
      rma<-ifelse(ord==0,upp,min(upp,sqrt(nn)))
      optimize(rhogam,lower=10^-10,upper=rma,gam0=gam0, r00=r,ord=ord,nn=nn)$objective
       }

getrgam<-function(gam0,ord,nn=0){
      if(gam0==0) getr0(ord,nn)
      else {rma<-ifelse(ord==0,upp,min(upp,sqrt(nn)))
            erg<-optimize(rhogamo,lower=10^-10,upper=rma,gam0=gam0,maximum=T,ord=ord,nn=nn)
            rho<-erg$objective; r0<-erg$maximum;
            if(ord==0) cr0<-c0opt(r0)
            if(ord==1) cr0<-c1opt(nn,r0)
            if(ord==2) cr0<-c2opt(nn,r0)
            c(rho=rho,r0=r0,cr0=cr0)
             }
}

upp<-100

####Examples
getrgam(gam0=0,ord=0) 
getrgam(gam0=0,ord=1,nn=10) 
getrgam(gam0=0,ord=2,nn=10) 
getrgam(gam0=2,ord=0) 
getrgam(gam0=2,ord=1,nn=10) 
getrgam(gam0=2,ord=2,nn=10) 
getrgam(gam0=3,ord=0) 
getrgam(gam0=3,ord=1,nn=10) 
getrgam(gam0=3,ord=2,nn=10) 


######################################################################################
#Table 1
# uses getoptpsi from FYZ.R
#      MSEn from MSEexact.R
######################################################################################

allMSEs<-function(r,n,AUS=F){ ## r: radius, n: sample size, AUS: mode of output
   if (n<50)
       METH<-cdfAlgC
   else 
       METH<-cdfAlgD     # n<30 => Alog C else Algo D
   #f-o
   c0<-c0opt(r)
   print("f-o")
   m0<-MSEn(cc=c0,n=n,r=r,fun=METH)$value  ## exact MSE of f-o-o M-est
   #s-o
   c1<-c1opt(n,r)
   print("s-o")
   m1<-MSEn(cc=c1,n=n,r=r,fun=METH)$value  ## exact MSE of s-o-o M-est
   #t-o
   c2<-c2opt(n,r)
   print("t-o")
   m2<-MSEn(cc=c2,n=n,r=r,fun=METH)$value  ## exact MSE of t-o-o M-est
   #FYZ
   tr<-try(getoptpsi(del=10^(-4),eps=r/sqrt(n),n=n),TRUE)  ## with try if optimization fails
   if(inherits(tr,"try-error"))  ## if getoptpsi fails -> NAs
      {cF<-NA
       mF<-NA}
   else
      {cF<-tr$c; 
       print("FYZ")
       mF<-MSEn(cc=cF,n=n,r=r,fun=METH)$value  ## exact MSE of FYZ M-est
       }
   #ex-opt
   print("ex")
   MSEex<-function(cx){MSEn(cc=cx,n=n,r=r,fun=METH)$value}      
   erg<-optimize(MSEex,interval=c(0.001,4))
   cE<-erg$minimum
   mE<-erg$objective
   ## relative performances
   re0<-(m0/mE-1)*100
   re1<-(m1/mE-1)*100
   re2<-(m2/mE-1)*100
   reF<-(mF/mE-1)*100  
   cat("\n")
   print(c(c0=c0,m0=m0,re0=re0,c1=c1,m1=m1,re1=re1,c2=c2,m2=m2,re2=re2,cF=cF,mF=mF,reF=reF,cE=cE,mE=mE),3)))
   if(AUS==T)
       list(c0=c0,m0=m0,re0=re0,c1=c1,m1=m1,re1=re1,c2=c2,m2=m2,re2=re2,cF=cF,mF=mF,reF=reF,cE=cE,mE=mE)
   else   c(c0=c0,m0=m0,re0=re0,c1=c1,m1=m1,re1=re1,c2=c2,m2=m2,re2=re2,cF=cF,mF=mF,reF=reF,cE=cE,mE=mE)
}
allMSEs(0.1,5)


#######################################################################################
# Cniper points
#######################################################################################

cniper<-function(n,r){
   c1<-c1opt(n,r)
    A<-Anorm(c1)
    V<-Vnorm(c1)
    b<-A*c1
    M<-asMSE0.norm(c1,r)
    cnip0<-sqrt(M-1)/r  ## f-o cniper-point
    cnip1<-sqrt(abs(M-1+(M+b^2*(r^2+1)+1)*r/sqrt(n))/(r^2*(1-1/n)+r/sqrt(n))) ## s-o cniper-point
    c(c1,cnip0,cnip1)
}

testrisks<-function(n){
    ord<-ifelse(n<10^12,1,0)
    r0<-getr0(ord,nn=n)[2]
    cnip0<-cniper(n,r0)[2]
#
    cnip<-cniper(n,r0)[3]
    pr<-1-pnorm(cnip)
    qn<-pr+r0/sqrt(n)*(1-pr)
     
   err.ii<-function(al,pr0=pr,qn0=qn,nn=n)
       { if(n<10^12)  ## determines error.II of Niveau alpha=a1 - NP Test
           { # bin-test niveau 0.05
             n0<-qbinom(1-al,size=nn,prob=pr0)
             n1<-n0-1
             pp0<-1-pbinom(n0,size=nn,prob=pr0)
             pp1<-1-pbinom(n1,size=nn,prob=pr0)
             gama<-(al-pp0)/(pp1-pp0) # randomization
             p<-1-pbinom(n0,size=nn,prob=qn0)+gama*dbinom(n0,size=nn,prob=qn0)
           }
         else
           { d<-qnorm(1-al)
             dn<-d-r0*sqrt((1-pr0)/pr0)
             p<-1-pnorm(dn)
           }
       return(1-p)
       } 

   error.ii<-err.ii(0.05)
   minmaxr<-function(a){a-err.ii(a)} 
   eps.n<-uniroot(minmaxr,low=0.00001,up=0.5)$root
   c(r=r0,cnip=cnip,cnip0=cnip0,error.ii=error.ii,eps.n=eps.n)
}