### Calculation of exact variance and MSE of the median               
# in the Gaussian location model                                      
#                                                                     
# Author : Dr. P. Ruckdeschel                                         
# Date   : 05.10.05                                                   
######################################################################

######################################################################
## in the ideal model:
######################################################################

#Attention:
# n in varmed, varmedmid, msemed, msemedmid stands for (n[-1])/2

## odd sample size

varmed<-function(n,bis=10){
   fdichte<-function(t,n){
     #density of the median
     d<-dnorm(t)
     p<-pnorm(t)
     t^2*(2*n+1)^2*exp(lgamma(2*n+1)-2*lgamma(n+1)+n*log(p*(1-p)))*d
     }
   2*integrate(fdichte,0,bis*log(2*n+1)/sqrt(2*n+1),n=n)$value
}

## even sample size --- midpoint estimator

varmedmid<-function(n,bis=10){
    # inner integrand -> gives H_0(u) in about (5.60)
    fidichte<-function(s,t,n){
       d1<-dnorm(s)
       d2<-dnorm(2*t-s)
       p1<-1-pnorm(s)
       p2<-pnorm(2*t-s)
       J0<-lgamma(2*n)-lgamma(n+1)-lgamma(n)
       exp(J0+(n-1)*log(p1*p2))*d1*4*n^2*d2
       }
    # outer integrand 
    fdichte<-function(t,n,bis0=bis){
       up<-ifelse(n<1000,15*n,10*n)
       L2f<-function(t)
           integrate(f=fidichte,lower=t,upper=t+log(up+1)/sqrt(n+1),
                     n=n,t=t,subdivisions=5000,rel.tol=10^(-10))$value
       # calculates H_0(t)
       L2<-sapply(t,L2f)
       t^2*(2*n)*L2
       }

     if(n>1)
        #integrate left an right part separately 
       {I1<-integrate(fdichte,-bis*log(2*n+1)^2/sqrt(2*n+1),0,
                      n=n,subdivisions=5000,rel.tol=10^(-10))$value
        I2<-integrate(fdichte,0, bis*log(2*n+1)^2/sqrt(2*n+1),
                      n=n,subdivisions=5000,rel.tol=10^(-10))$value
        return(I1+I2)}
     else
       {return(1)}
    }

######################################################################
## in the contaminated model:
######################################################################

## odd sample size

msemedodd<-function(n,r=0.1,bis=10){

 # integrand t^2 * g_{n,k,k} from (5.8)
    fkdichte<-function(t,n,k1){
       d<-dnorm(t)
       p<-pnorm(t)
       return(t^2*(2*n+1)*(2*n+1-k1)*exp(lgamma(2*n+1-k1)-
              lgamma(n+1-k1)-lgamma(n+1)+n*log(p)+(n-k1)*log(1-p))*d)
                              }
 # gives n MSE(med|K=k)
    varquant<-function(k0,n0,bis0=10){
            integrate(fkdichte,-bis0*2/sqrt(2*n0),bis0*2/sqrt(2*n0),n=n0,k1=k0)$value}

 # calculating censored binomial probabilities 
    vvec<-sapply(0:n,varquant,n0=n,bis0=bis)
    pvec<-dbinom(0:n,size=2*n+1,prob=r/sqrt(2*n+1))/
          pbinom(q=n,size=2*n+1,prob=r/sqrt(2*n+1))
    evec<-sum(pvec*(0:n))
 # gives n MSE(med)
    sum(vvec*pvec)
}

## even sample size --- midpoint estimator


msemedmid<-function(n,r=0.1,bis=10){
 # inner integrand 
    fidichte<-function(s,ui,ni,ki){
       ds<-dnorm(s)
       ps<-pnorm(s)
       (ui+s)^2*ps^(ni-1-ki)*ds
       }
 # inner integral (H(u)) to be evaluated pointwise
    Huf<-function(u,n,k,up=3){
         integrate(f=fidichte,lower=u-2*log(up)/sqrt(n+1),
                   upper=u,ni=n,ui=u,ki=k,subdivisions=5000,
                   rel.tol=10^(-10))$value}
 # inner integral (H(u)) for vector-wise evaluation
    Hu<-function(u,n,k,bis0=40){
        up0<-ifelse(n<1000,15*n,10*n)
        sapply(u,Huf,n=n,k=k,up=up0)
        }
 # ounter integrand 
    fkdichtee<-function(u,n,k1,aus="nein"){
        d<-dnorm(u)
        p<-pnorm(u)
        Hu(u=u,n=n,k=k1)*d*exp((n-1)*log(1-p)+
        lgamma(2*n-k1)-lgamma(n+1-k1)-lgamma(n))*n/2*(2*n-k1)*(n-k1)
        }
 # gives n MSE(med|K=k)
    varquante<-function(k0,n0,bis0=10){
        sds=ifelse(n0<100,5000,4000);tol=ifelse(n0<100,10^(-10),10^(-8))
        I1<-integrate(fkdichtee,low=0,up=bis0*log(2*n0-1.5*k0)/sqrt(2*n0),
                      n=n0,k1=k0,subdivisions=sds,rel.tol=tol)$value
        I2<-integrate(fkdichtee,low=-bis0*log(2*n0-1.5*k0)/sqrt(2*n0),up=0,
                      n=n0,k1=k0,subdivisions=sds,rel.tol=tol)$value
        I1+I2;
        }
 # calculating censored binomial probabilities 
     vvec<-sapply(0:n-1,varquante,n0=n,bis0=bis)
     pvec<-dbinom(0:n-1,size=2*n,prob=r/sqrt(2*n))/pbinom(q=n-1,size=2*n,prob=r/sqrt(2*n))
     evec<-sum(pvec*(0:n-1))
 # gives n MSE(med)
     sum(vvec*pvec)
}


################################################################
#envelope  function for all cases
################################################################
## returns numerically exact value ("tat") and asymptotics (0-2)

msemed<-function(n,r){
if(n<=2)
{ c("as0"=pi/2*(1+r^2),
    "as1"=pi/2*(1+r^2)*(1+2*r/sqrt(n)),
    "as2"=pi/2*((1+r^2)*(1+2*r/sqrt(n))+1/n*(3*r^4+3*r^2-3+pi*(r^4+6*r^2+3)/6)),
    "tat"=1)}
else
 {          # the integration bound is a bit tricky....
  if (n==3) bisa=5
  if (n==4) bisa=10
  if (n==5) bisa=8
  if (n>5)  bisa=10
  if (r==0)
      {
       if(n%%2==0)
          { c("as0"=pi/2,
              "as1"=pi/2,
              "as2"=pi/2*(1+1/n*(-3+pi/2)),
              "tat"=varmedmid((n/2),bis=bisa))          
          }
       else
          { c("as0"=pi/2,
              "as1"=pi/2,
              "as2"=pi/2*(1+1/n*(-2+pi/2)),
              "tat"=varmed((n-1)/2,bis=bisa))
          }
      }
  else{
       if(n%%2==0)
          { if(n<=500)
              c("as0"=pi/2*(1+r^2),
                "as1"=pi/2*(1+r^2)*(1+2*r/sqrt(n)),
                "as2"=pi/2*((1+r^2)*(1+2*r/sqrt(n))+1/n*(3*r^4+3*r^2-3+pi*(r^4+6*r^2+3)/6)),
                "tat"=msemedmid((n/2),r,bis=bisa))
            else ## exact value only if n<=500 !
              c("as0"=pi/2*(1+r^2),"as1"=pi/2*(1+r^2)*(1+2*r/sqrt(n)),
                "as2"=pi/2*((1+r^2)*(1+2*r/sqrt(n))+1/n*(3*r^4+3*r^2-3+pi*(r^4+6*r^2+3)/6)),
                "tat"=NA)
          }
       else
          { c("as0"=pi/2*(1+r^2),
             "as1"=pi/2*(1+r^2)*(1+2*r/sqrt(n)),
             "as2"=pi/2*((1+r^2)*(1+2*r/sqrt(n))+1/n*(3*r^4+3*r^2-2+pi*(r^4+6*r^2+3)/6)),
             "tat"=msemedodd((n-1)/2,r,bis=bisa))
          }
      }
  }
}      

relerror<-function(n,r){
   A<-msemed(n,r)
   100*c("as0"=A[1]/A[4]-1,
         "as1"=A[2]/A[4]-1,
         "as2"=A[3]/A[4]-1,
         "eps"=r/sqrt(n)/100)
     }