################################################################################ ## MSE -optimal IC according to FYZ(01) # (c) P. Ruckdeschel 24-02-04 ################################################################################ la<-function(x,mu){tanh(x*mu)} #lambda in paper (4.2) de<-function(x,mu){x-mu*tanh(x*mu)} #delta in paper lac<-function(x,mu,cc){x0<-pmin(cc,pmax(-cc,x)); la(x0,mu)} #clipped versions (4.1) dec<-function(x,mu,cc){x0<-pmin(cc,pmax(-cc,x)); de(x0,mu)} De<-function(x,mu){dnorm(mu+x)+dnorm(mu-x)} #(4.4) psi<-function(x,mu,a,b,cc){ ### optimal type of IC x0<-pmin(cc,pmax(-cc,x)) a*la(x0,mu)+b*de(x0,mu) } getmuo<-function(del=10^(-6),eps,muu){ ### to given eps determines \bar mu (4.8) Im<-function(mu,eps,del){ Imi<-function(x,mu){la(x,mu)^2*De(x,mu)} integrate(Imi,0,Inf,rel.tol=del,mu=mu)$value-eps/(1-eps) } te<-uniroot(Im,c(muu,10),tol=del,del=del,eps=eps)$root min(te,1) } c1mu<-function(mu,del=10^(-6),eps){ ###to given mu,eps determines clipping height c1 in (A.40) Ic1m<-function(cc,mu=mu,del=del,eps=eps){ Ic1i<-function(x,cc,mu){lac(x=x,mu=mu,cc=cc)*la(x=x,mu=mu)*De(x=x,mu=mu)} integrate(Ic1i,0,Inf,cc=cc,mu=mu,rel.tol=del)$value-eps/(1-eps)*lac(x=cc,mu=mu,cc=cc) } uniroot(Ic1m,c(10^(-9),10000),tol=del,del=del,eps=eps,mu=mu)$root } c2mu<-function(mu,del=10^(-6),eps){ ###to given mu,eps determines clipping height c2 in (A.41) Ic2m<-function(cc,mu,del=10^(-6),eps){ Ic2i<-function(x,cc,mu){dec(x,mu,cc)*la(x,mu)*De(x,mu)} integrate(Ic2i,0,Inf,cc=cc,mu=mu,rel.tol=del)$value-eps/(1-eps)*dec(cc,mu,cc) } uniroot(Ic2m,c(10^(-9),10000),tol=del,del=del,eps=eps,mu=mu)$root } K1mu<-function(mu,del=10^(-6),eps){ ###to given mu,eps determines K1 in (A.40) c1<-c1mu(mu,del,eps=eps) IK1i<-function(x,cc,mu){lac(x,mu,cc)*de(x,mu)*De(x,mu)} IK1m<-integrate(IK1i,0,Inf,cc=c1,mu=mu,rel.tol=del)$value la(c1,mu)/IK1m } K2mu<-function(mu,del=10^(-6),eps){ ###to given mu,eps determines K2 in (A.41) c2<-c2mu(mu,del,eps=eps) IK2i<-function(x,cc,mu){dec(x,mu,cc)*de(x,mu)*De(x,mu)} IK2m<-integrate(IK2i,0,Inf,cc=c2,mu=mu,rel.tol=del)$value de(c2,mu)/IK2m } std1mu<-function(c1,mu,del=10^(-6),eps){ ###to given mu,c1,eps determines sigma1 for getvar1 IK1i<-function(x,cc,mu){lac(x,mu,cc)*de(x,mu)*De(x,mu)} IK1m<-integrate(IK1i,0,Inf,cc=c1,mu=mu,rel.tol=del)$value IK1m } std2mu<-function(c2,mu,del=10^(-6),eps){ ###to given mu,c2,eps determines sigma2 for getvar2 IK2i<-function(x,cc,mu){dec(x,mu,cc)*de(x,mu)*De(x,mu)} IK2m<-integrate(IK2i,0,Inf,cc=c2,mu=mu,rel.tol=del)$value IK2m } getab<-function(cc,mu,del=10^(-6),eps){ ###to given mu,cc,eps determines a,b according to display before (A.45) I1i<-function(x,mu,cc){dec(x,mu,cc)*de(x,mu)*De(x,mu)} I2i<-function(x,mu,cc){lac(x,mu,cc)*de(x,mu)*De(x,mu)} I3i<-function(x,mu,cc){dec(x,mu,cc)*la(x,mu)*De(x,mu)} I4i<-function(x,mu,cc){lac(x,mu,cc)*la(x,mu)*De(x,mu)} I1<-integrate(I1i,0,Inf,cc=cc,mu=mu,rel.tol=del)$value I2<-integrate(I2i,0,Inf,cc=cc,mu=mu,rel.tol=del)$value I3<-integrate(I3i,0,Inf,cc=cc,mu=mu,rel.tol=del)$value I4<-integrate(I4i,0,Inf,cc=cc,mu=mu,rel.tol=del)$value dc<-de(cc,mu) lc<-la(cc,mu) A11<-I1 A12<-I2 A21<-I3-dc*eps/(1-eps) A22<-I4-lc*eps/(1-eps) a<- -A21/(A11*A22-A12*A21) b<- A22/(A11*A22-A12*A21) #print(c(a=a,b=b,c=cc,mu=mu)) return(list(a=a,b=b)) } getvar<-function(cc,mu,del=10^(-6),eps,env){ ###to given mu,cc,eps determines variance of psi acc. to (4.6) erg<-getab(cc,mu,del,eps) a<-erg[[1]] b<-erg[[2]] assign("psi0",function(x){psi(x,mu,a,b,cc)},pos=env) ###writes psi0, a, b to the calling environment assign("a",a,pos=env) assign("b",b,pos=env) dc<-de(cc,mu) lc<-la(cc,mu) Ii<-function(x,mu,aa,bb,cc){psi(x,mu=mu,a=aa,b=bb,cc=cc)^2*De(x,mu)} I<-integrate(Ii,0,Inf,mu=mu,aa=a,bb=b,cc=cc,rel.tol=del)$value return(((1-eps)*I+eps*(a*lc+b*dc)^2)/(1-eps)^2) } getvar1<-function(mu,del=10^(-6),eps,env){ ###to given mu,eps determines variance of lac acc. to (4.6) c1<-c1mu(mu=mu,del=del,eps=eps) assign("cc1",c1,pos=env) ###writes psi1, cc1 to the calling environment lc<-la(c1,mu) Ii<-function(x,mu,cc){lac(x,mu=mu,cc=cc)^2*De(x,mu)} I<-integrate(Ii,0,Inf,mu=mu,cc=c1,rel.tol=del)$value std<-std1mu(c1,mu,del,eps) assign("psi1",function(x){lac(x,mu=mu,cc=c1)/std},pos=env) return(((1-eps)*I+eps*lc^2)/(1-eps)^2/std^2) } getvar2<-function(mu,del=10^(-6),eps,env){ ###to given mu,eps determines variance of dec acc. to (4.6) c2<-c2mu(mu,del,eps) assign("cc2",c2,pos=env) ###writes psi2, cc2 to the calling environment dc<-de(c2,mu) Ii<-function(x,mu,cc){dec(x,mu=mu,cc=cc)^2*De(x,mu)} I<-integrate(Ii,0,Inf,mu=mu,cc=c2,rel.tol=del)$value std<-std2mu(c2,mu,del,eps) assign("psi2",function(x){dec(x,mu=mu,cc=c2)/std},pos=env) return(((1-eps)*I+eps*dc^2)/(1-eps)^2/std^2) } getmse<-function(cc,mu,del=10^(-6),eps,n,sit,env,varmse=F){ #to given mu,(cc,)eps determines one of the competing MSEs if(sit==0) v<-getvar(cc,mu,del,eps,env)/n if(sit==1) v<-getvar1(mu,del,eps,env)/n if(sit==2) v<-getvar2(mu,del,eps,env)/n if(varmse==T) ##returns var if varmse==T else MSE return(v) else return(mu^2+v) } getoptpsi<-function(del=10^(-6),eps,n,varmse=FALSE,aus=FALSE,graf=FALSE){ muu<-qnorm(1/2/(1-eps)) # \underline \mu acc. to (A.34) muo<-getmuo(del,eps,muu) cc0<-0 cc1<-0 cc2<-0 envi<-environment() MSEc<-function(mu,del,eps,n,env,aus){ ## optimizes cc to given mu in situation (1) ccu0<-c1mu(mu,del,eps) cco0<-c2mu(mu,del,eps) ergc0<-optimize(getmse,interval=c(min(ccu0,cco0),max(ccu0,cco0)),maxim=F,tol=del,mu=mu,del=del,eps=eps,n=n,env=envi,sit=0,varmse=varmse) cc0<-ergc0$minimum mse0<-ergc0$objective assign("cc0",cc0,pos=env) if(aus==T) print(c(muu=muu,muo=muo,mu=mu,cu=min(ccu0,cco0),co=max(ccu0,cco0),c0=cc0,MSE0=mse0*n)) return(mse0) } erg0<-optimize(MSEc,interval=c(muu,muo),maxim=F,tol=del,del=del,eps=eps,n=n,env=envi,aus=aus) ## optimizes mu in situation (0) # cc0 in this environment is written from the function MSEc mu0<-erg0$minimum; mse0<-erg0$objective if(aus==T) print(c(a=a,b=b,mu0=mu0,c0=cc0,MSE0=mse0*n)) cc<-cc0 mu<-mu0 mse<-mse0 psio<-psi0 # erg1<-optimize(getmse,interval=c(muu,muo),maxim=F,tol=del,del=del,eps=eps,n=n,env=envi,sit=1,cc=NA,varmse=varmse) ## optimizes mu in situation (1) # cc1 in this environment is written from the function MSEc mu1<-erg1$minimum; mse1<-erg1$objective if(aus==T) print(c(a=1,b=0,mu1=mu1,c1=cc1,MSE1=mse1*n)) if(mse1