#############M-Estimators #MAPLE files for #[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 # # ######################################### # conventions: ######################################### # t ^= $t^\natural$ from after (A.27) # l2 ^= $l_2$ from (D) # l3 ^= $l_3$ from (D) # v0 ^= $v_0$ from (D) # v1 ^= $\tilde v_1$ from (D) # v2 ^= $\tilde v_2$ from (D) # r0 ^= $\rho_0$ from (D) # r1 ^= $\rho_1$ from (D) # k0 ^= $\kappa_0$ from (D) # w ^= $\sqrt(\bar n)$ from (A.6) # S ^= $\tilde s(u)$ from after (A.27) # SZ ^= numerator of S # SN ^= denominator of S # S1 ^= $s'(u)$ # SZ1 ^= numerator of S1 # SN1 ^= denominator of S1 # G10 ^= $G^{(1)}_{n,u/\sqrt{\bar n}}$ from (A.32) # G20 ^= $G^{(2)}_{n,u/\sqrt{\bar n}}$ from after (A.32) # K ^= $k^\natural$ (A.41) # k ^= $\tilde k$ (A.41) # W ^= $\sqrt(n)$ # i ^= order of approximation : # 3 ^= up to o(n^{-1}) (t-o) (maximum!) # 2 ^= up to o(n^{-1/2}) (s-o) # 1 ^= up to o(n^{0}) (f-o) # pm ^= positive or negative contamination: # 1 ^= positive # -1 ^= negative #[Ru1] ##calculation of asymptotic expansion of S asS:=proc(i, withprint) local SZ, SN, S; SZ:=u-t-l2*u^2/w/2-l3*u^3/w^2/6; SN:=v0*(1+v1*u/w+v2*u^2/2/w^2); S:=collect(subs(w=1/w,convert(subs(w=1/w,asympt(SZ/SN,w,i)),polynom)),w); if withprint=1 then print (S); fi; RETURN(eval(S)) end; #[Ru1] ##calculation of asymptotic expansion of S1 asS1:=proc(i, withprint) local SZ1, SN1, S1, SZ, SN; SZ:=u-t-l2*u^2/w/2-l3*u^3/w^2/6; SN:=v0*(1+v1*u/w+v2*u^2/2/w^2); SZ1:=diff(SZ,u); SN1:=diff(SN,u); S1:=collect(subs(w=1/w,convert(subs(w=1/w,asympt(SZ1/SN-SZ*SN1/SN^2,w,i)),polynom)),w); if withprint=1 then print (S1); fi; RETURN(eval(S1)) end; #[Ru1] ##calculation of asymptotic expansion of $\tilde g_n$ from (A.33) asg:=proc(i, withprint) local G10, G20, S1, S, SS, gn0, gn, gn1, gn2, Ph, P0,P1, P2; ## to ease computations: # s^= asy up to 1/w^2, # ss^= asy up to 1/w, ## G10:=1+r0/(6*w)*(s^3-3*s)+r1*u/(6*w^2)*(ss^3-3*ss)+k0*(ss^4-6*ss^2+3)/(w^2*24)+r0^2*(ss^6-15*ss^4+45*ss^2-15)/(w^2*72); G20:=r1*(1-ss^2)/w^2/6; S:=asS(i,0); SS:=asS(i-1,0); S1:=asS1(i,0); # $\tilde g_n(u)$ from (8.44) gn0:=collect(subs(w=1/w,convert(subs(w=1/w,asympt(collect(v0*(G10*S1+G20),w),w,i)),polynom)),w); # substitute $s,ss=\tilde s(u)$ gn:=collect(collect(subs(w=1/w,convert(subs(w=1/w,asympt(collect(subs(s=S,ss=SS,gn0),w),w,i)),polynom)),v0),w); if withprint=1 then gn1:= collect(collect(subs(w=1/w,convert(subs(w=1/w,asympt(collect(gn,w),w,1)),polynom)),v0),w); gn2:= collect(collect(subs(w=1/w,convert(subs(w=1/w,asympt(collect(gn,w),w,2)),polynom)),v0),w); P2:=collect(collect(collect(collect(gn-gn2,u),r0),v0),w); P1:=collect(collect(collect(collect(gn2-gn1,u),r0),v0),w); print ("P1"); print(P1); print ("P2"); print(P2); fi; RETURN(eval(gn)) end; asg(3,1); #[Ru1] ##ds^= $\tilde s- s_1$ in (8.53) ##dfac ^= $\varphi(\tilde s)/ \varphi(\tilde s_1)$ in (A.35) #[Ru1] ##calculation of asymptotic expansion of dfac asdfac:=proc(i, withprint) local s1, ds, dss, Ds, Dss, dfac, dfac0, dfac1; ## to ease computations: # Ds, ds ^= asy up to 1/w^2, # Dss, dss^= asy up to 1/w, ## s1:=(u-t)/v0; ds:=asS(i,0)-s1; dss:=asS(i-1,0)-s1; # dfac0:=1-s*Ds+(s^2-1)*Dss^2/2; dfac1:= subs(Ds=ds,Dss=dss,dfac0); dfac:=collect(subs(w=1/w,convert(subs(w=1/w,asympt(collect(dfac1,w),w,i)),polynom)),w); if withprint=1 then print (dfac); fi; RETURN(eval(dfac)) end; #[Ru1] ##calculation of asymptotic expansion of $g_n$ from before (A.36) asgns:=proc(i, withprint) local u1, dfacs, gns0, g, g2, g1, P1, P2; u1:=v0*s+t; dfacs:=collect(collect(collect(subs(u=u1,asdfac(i,0)),v0),s),w); gns0:=collect(collect(collect(subs(u=u1,asg(i,0)),v0),s),w); # g:=collect(collect(collect(subs(w=1/w,convert(subs(w=1/w,asympt(collect(collect(collect(dfacs*gns0,v0),s),w),w,i)),polynom)),t),v0),w); if withprint=1 then print (g); g1:= collect(collect(collect(subs(w=1/w,convert(subs(w=1/w,asympt(collect(g,w),w,1)),polynom)),t),v0),w); g2:= collect(collect(collect(subs(w=1/w,convert(subs(w=1/w,asympt(collect(g,w),w,2)),polynom)),t),v0),w); P2:=g-g2; P1:=g2-g1; print ("tilde P1"); print(P1); print ("tilde P2"); print(P2); fi; RETURN(eval(g)) end; asgns(3,1); #[Ru1] ##calculation of asymptotic expansion of $h_n$ from (A.37)/(A.38) ashn:=proc(i, pot, withprint) local u1, gn, hn, hn0, h1, h2, Q1,Q2; u1:=v0*s+t; gn:=asgns(i,0); hn0:=gn*u1^pot; # hn:=collect(collect(collect(collect(subs(w=1/w,convert(subs(w=1/w,asympt(collect(collect(collect(hn0,v0),s),w),w,i)),polynom)),v0),t),s),w); if withprint=1 then print (hn); h1:= collect(collect(collect(collect(subs(w=1/w,convert(subs(w=1/w,asympt(collect(hn,w),w,1)),polynom)),v0),t),s),w); h2:= collect(collect(collect(collect(subs(w=1/w,convert(subs(w=1/w,asympt(collect(hn,w),w,2)),polynom)),v0),t),s),w); Q2:=hn-h2; Q1:=h2-h1; print ("Q1"); print(Q1); print ("Q2"); print(Q2); fi; RETURN(eval(hn)) end; #[Ru1] ###least favorable contamination: (A.39) #[Ru1] #sort by t HN:=collect(collect(collect(collect(ashn(3,2,0),v0),s),t),w); #differentiate: HND:=diff(HN,t); HND2:=collect(collect(collect(diff(HND,t),t),s),w); HND2s:=collect(collect(collect(subs(r0=0,v1=0,l2=0,HND2),t),s),w); # K ^= $K^natural$ (A.41) ##=> t= +-Kb #[Ru1] ##calculation of asymptotic expansion of $\tilde h_n$ from (A.42) ash:=proc(i, pot, pm, withprint) local hn0,hn; # hn0:=collect(collect(collect(collect(ashn(i,pot,0),v0),s),t),w); hn:=collect(collect(collect(collect(subs(w=1/w,convert(subs(w=1/w,asympt(collect(collect(collect(subs(t=pm*K*b,hn0),v0),s),w),w,i)),polynom)),v0),t),s),w); if withprint=1 then print (hn); fi; RETURN(eval(hn)) end; #[Ru1] ### integration w.r.t. s (A.45) intesout:=proc(expr,withprint) local tt1,tt2,tt3,tt4,tt,btt; # ### Moments of N(0,1) (A.44) M2:=1; M4:=3; M6:=15; M8:=105; # tt1:=subs(s^8=M8,collect(collect(expr,s),w)); tt2:=subs(s^6=M6,tt1); tt3:=subs(s^4=M4,tt2); tt4:=subs(s^2=M2,tt3); tt:=collect(collect(subs(s=0,tt4),k),w); if withprint=1 then print(tt); fi; RETURN(eval(tt)) end; ##[Ru1] asMSEK:=proc(i, pm, withprint) local hn0,hn, MSEK, M2, M4, M6, M8; # hn0:=ash(i,2,pm,0); MSEK:=collect(collect(collect(collect(collect(intesout(hn0,0),v1),v0),b),K),w);; if withprint=1 then print (MSEK) fi; RETURN(eval(MSEK)) end; ### collection of terms ## N-nat in terms of K-tilde asNn:=proc(i, withprint) local nnat0, nnat; ## (A.19) nnat0:=sqrt(1-k/W); nnat:=collect(collect(subs(W=1/W,convert(subs(W=1/W,asympt(nnat0,W,i)),polynom)),k),W); if withprint=1 then print (nnat); fi; RETURN(eval(nnat)) end; ## K-nat in terms of K-tilde asKn:=proc(i, withprint) local knat0, knat; ## (A.41) knat0:=k/asNn(i,0); knat:=collect(collect(subs(W=1/W,convert(subs(W=1/W,asympt(knat0,W,i)),polynom)),k),W); if withprint=1 then print (knat) fi; RETURN(eval(knat)) end; ## (A.49) asMSEk:=proc(i, pm, withprint) local fact, nnat, k0, MSEK, MSEk, MSEk0; # fact:=1/asNn(i,0)^2; MSEK:=asMSEK(i,pm,0); nnat:=asNn(i,0); k0:=asKn(i,0); MSEk0:=collect(collect(subs(K=k0,fact*subs(w=W*nnat,MSEK)),k),W); MSEk:=collect(collect(collect(collect(subs(W=1/W,convert(subs(W=1/W,asympt(MSEk0,W,i)),polynom)),v0),b),k),W); if withprint=1 then print (MSEk) fi; RETURN(eval(MSEk)) end; ##Moments of (W K)~Bin(W^2,r/W) Binmoment:=proc(i) local KKi,p,x,KKis; KKi:=x^i*binomial(n,x)*p^x*(1-p)^(n-x); KKis:=collect(subs(p=r/W,n=W^2,simplify(subs((-1/(-1+p))^n=1/(1-p)^n,simplify(sum(KKi,x=0..n))))/W^i),W); RETURN(eval(KKis)) end; Binmoment(4); for i to 4 do K[i] := Binmoment(i); od; intekout:=proc(expr,i) local expr0, j0,KK, expro; expr0:=collect(expr,W); for j0 to 4 do KK := Binmoment(5-j0); expr0:=subs(k^(5-j0)=KK,expr0); od; expro:=collect(subs(W=1/W,convert(subs(W=1/W,asympt(expr0,W,i)),polynom)),W); RETURN(eval(expro)) end; asMSE:=proc(i, pm, withprint) local MSEk, MSE0, MSE; # MSEk:=asMSEk(i,pm,0); MSE0:=intekout(MSEk,i); MSE:=collect(collect(collect(collect(MSE0,b),v0),r),W); if withprint=1 then print (MSE) fi; RETURN(eval(MSE)) end; ## expectation of S^i under Qn asESi:=proc(i,pot, pm, withprint) local fact, nnat, k0, hn0, ESiK, ESik, ESik0, ESi0, ESi; # hn0:=ash(i,pot,pm,0); ESiK:=collect(collect(collect(collect(collect(intesout(hn0,0),v1),v0),b),K),w);; # fact:=1/asNn(i,0)^pot; nnat:=asNn(i,0); k0:=asKn(i,0); ESik0:=collect(collect(subs(K=k0,fact*subs(w=W*nnat,ESiK)),k),W); ESik:=collect(collect(collect(collect(subs(W=1/W,convert(subs(W=1/W,asympt(ESik0,W,i)),polynom)),v0),b),k),W); ESi0:=intekout(ESik,i); ESi:=collect(collect(collect(collect(ESi0,b),v0),r),W); if withprint=1 then print (ESi) fi; RETURN(eval(ESi)) end; ### control MSE1:=asMSE(3,1,0);MSE2:=asESi(3,2,1,0); MSE1-MSE2; ##(A.49)-(A.53) M1:=asMSE(3,1,0); M2:=asMSE(3,-1,0); M3:=asMSE(3,0,0); ### decision if sup psi = -inf psi: dM:=collect(collect(collect(collect(W/r/b*(M1-M2)/2,v0),b),r),v1); dM1:=diff(dM,v1); dM0:=subs(v1=0,dM); V10:=collect(collect(collect(collect(subs(W=1/W,convert(subs(W=1/W,asympt(4/l2*dM0/dM1,W,3)),polynom)),v0),b),r),W)*l2/4; V1:=-l2/4*(b^2/v0^2*(r^2+3)*(1+r/W-2*r^2/W^2)+3*(1-b^2/v0^2)); #(b^2/v0^2*r^2+3)*(1+r/W); #control: MM1:=collect(collect(collect(collect(subs(v1=V1,M1),v0),b),r),W); MM2:=collect(collect(collect(collect(subs(v1=V1,M2),v0),b),r),W); collect(collect(collect(collect(collect(MM1-MM2,l2),v0),b),r),W); B1:=asESi(3,1,1,1); B2:=asESi(3,1,-1,1); C1:=collect(collect(collect(collect(subs(W=1/W,convert(subs(W=1/W,asympt(B1^2,W,3)),polynom)),b),v0),r),W); C2:=collect(collect(collect(collect(subs(W=1/W,convert(subs(W=1/W,asympt(B2^2,W,3)),polynom)),b),v0),r),W); V1:=collect(collect(collect(collect(M1-C1,b),v0),r),W); V2:=collect(collect(collect(collect(M2-C2,b),v0),r),W); ##################### # #[Ru2]-Thm. 3.1 # ##################### getSa:=proc(a,i, withprint) local Sa0, Nn, Sa; Sa0:=subs(u=a,asS(i,0)); Nn:=asNn(i,0); Sa:=collect(subs(W=1/W,convert(subs(W=1/W,asympt(subs(w=W*Nn,subs(t=t/Nn,Sa0)*Nn),W,i)),polynom)),w); if withprint=1 then print (Sa); fi; RETURN(eval(Sa)) end; getdSa:=proc(a,i, withprint) local Sa0, Sa1, dSa; Sa0:=getSa(a,1,0); Sa1:=getSa(a,i,0); dSa:=Sa1-Sa0; if withprint=1 then print (dSa); fi; RETURN(eval(dSa)) end; ## asds:=proc(a,i,pm, withprint) local ds, Ds, dfac, dfac0, dfac1; ds:=getdSa(a,i,0); # dfac0:=Phis+pm*phis*(rho0/W/6*(s^2-1)+Ds); dfac1:= subs(s=getSa(a,1,0),subs(Ds=ds,dfac0)); dfac:=collect(subs(W=1/W,convert(subs(W=1/W,asympt(collect(dfac1,W),W,i)),polynom)),W); if withprint=1 then print (dfac); fi; RETURN(eval(dfac)) end; asdsr:=proc(i,pm,withprint) local dsm, dsp, dsrp, dsrm, dsr,k0; k0:=asKn(i,0); dsm:=asds(-a1,i,-1,0); dsp:=asds(a2,i,+1,0); # dsrm:=collect(subs(W=1/W,convert(subs(W=1/W,asympt(collect(subs(k=r,subs(K=k0,subs(t=K*bm,dsm))),W),W,i)),polynom)),W); dsrp:=collect(subs(W=1/W,convert(subs(W=1/W,asympt(collect(subs(k=r,subs(K=k0,subs(t=-K*bp,dsp))),W),W,i)),polynom)),W); if withprint=1 then print (dsrp); print (dsrm); fi; if (pm=1) then dsr:=dsrp else dsr:=dsrm fi; RETURN(eval(dsr)) end; asdsr(2,1,1);