### 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) }