NMST 434, Exercise session VI: M-estimators
March 26, 2020

SetOptions[Plot,BaseStyle->FontSize->16];

Example: (Subscript[S, n]) (sample standard deviation) as an M-estimator

Clear[M];
$Assumptions=τ>0 && σ2>0;
ψ[μ_,σ2_,τ_]={x-μ,(((x-μ))^2)-σ2,√σ2-τ}; (* esimating function ψ for the Z-estimator *)
D[ψ[μ,σ2,τ],{{μ,σ2,τ}}]/.x->μ (* matrix Γ *)
Ai=Inverse[%]/.σ->√σ2;
B=((Outer[Times,ψ[μ,σ2,τ],ψ[μ,σ2,τ]]//Expand)/.Table[(x^(4-j))->M[4-j],{j,0,4}]) (* matrix Σ with M[j] the j-th moment of X *)
Ai.B.Transpose[Ai]/.{M[1]->μ,M[2]-> σ2+(μ^2),σ2->(τ^2)}//Simplify//MatrixForm (* asymptotic variance matrix *)

Direct use of a Δ-theorem (for comparison, must give the same result)

Clear[M]
g[a_,b_]={a,b-(a^2),√(b-(a^2))}; (* g function to use in the Δ-theorem *)
(J=(D[g[a,b],{{a,b}}]/.{a->M[1],b->M[2]}))//MatrixForm (* Jacobian matrix of function g *)
Σ=Table[M[i+j]-M[i]M[j],{i,1,2},{j,1,2}]; (* variance matrix of (X,X^2) *)
J.Σ.Transpose[J]//Simplify//MatrixForm
Ai.B.Transpose[Ai]/.{μ->M[1],σ2->M[2]-((M[1])^2),τ->√(M[2]-((M[1])^2))}//Simplify//MatrixForm
(* the two as. variance matrices are indeed the same *)

Special case of normal distribution

M[k_]=Moment[NormalDistribution[μ,σ],k];
J.Σ.Transpose[J]//Simplify//MatrixForm
Ai.B.Transpose[Ai]/.{μ->M[1],σ2->M[2]-((M[1])^2),τ->√(M[2]-((M[1])^2))}//Simplify//MatrixForm

Example: Sample correlation coefficient as an M-estimator

For a general, properly integrable distribution. We find the joint distribution of the sample means of X and Y, sample variances of X and Y, sample average of XY, and the sample correlation coefficient.
For the joint distribution, WLOG centered distributions with unit variances does not apply (compare with the example Session1.nb). We must compute in full generality.

Clear[M];
$Assumptions=Element[μx,Reals]&&Element[μy,Reals]&&Element[μxy,Reals]&&σx2>0&&σy2>0&&-1<ρ<1;
ψ[μx_,μy_,μxy_,σx2_,σy2_,ρ_]={x-μx,y-μy,x y - μxy,(((x-μx))^2)-σx2,(((y-μy))^2)-σy2,((μxy - μx μy)/(√(σx2 σy2)))-ρ};
(A0 = D[ψ[μx,μy,μxy,σx2,σy2,ρ],{{μx,μy,μxy,σx2,σy2,ρ}}]);
(A = A0/.{x->M[1,0],y-> M[0,1]});
(A=A/.μxy->M[1,1])//MatrixForm
(B0 = Outer[Times,ψ[μx,μy,μxy,σx2,σy2,ρ],ψ[μx,μy,μxy,σx2,σy2,ρ]]);
(B=(B0//Expand)/.Flatten[Table[(x^(4-j))(y^(4-k))->M[4-j,4-k],{j,0,4},{k,0,4}]]);
(B=B/.μxy->M[1,1]//Simplify)//MatrixForm

Asymptotic variance matrix for general distributions.

(Var=Inverse[A].B.Transpose[Inverse[A]]//Simplify)//MatrixForm
Var/.{μx->M[1,0],μy->M[0,1],σx2->M[2,0]-((M[1,0])^2),σy2->M[0,2]-((M[0,1])^2)}//Simplify//MatrixForm

Asymptotic variance of the sample correlation coefficient.

Clear[M]
Short[(Var[[6,6]]/.{μx->M[1,0],μy->M[0,1],σx2->M[2,0]-((M[1,0])^2),σy2->M[0,2]-((M[0,1])^2)})//Simplify]
Var[[6,6]]/.{μx->M[1,0],μy->M[0,1],σx2->M[2,0]-((M[1,0])^2),σy2->M[0,2]-((M[0,1])^2),ρ->0}//Simplify

Special case of the bivariate normal distribution.

dist = MultinormalDistribution[{μx,μy},{{σx2,ρ √(σx2 σy2)},{ρ √(σx2 σy2),σy2}}];
M[i_,j_]=Moment[dist,{i,j}];
Var//Simplify//MatrixForm
Var/.{μx->0,μy->0,σx2->1,σy2->1}//Simplify//MatrixForm

Direct use of the Δ-theorem (only for the asymptotic distribution of (OverHat[ρ]), see Session1.nb)

Clear[μ];
g[a_,b_,c_,d_,e_]=(((e-a b))/((((c-(a^2))(d-(b^2))))^(1/2)));
Dg = D[g[a,b,c,d,e],{{a,b,c,d,e}}]/.{a->0,b->0,c->1,d->1,e->ρ}//Simplify;
MatrixForm[Dg](* gradient of g at E(a,b,c,d,e) *)
Σ={{1,ρ,μ[3,0],μ[1,2],μ[2,1]},{ρ,1,μ[2,1],μ[0,3],μ[1,2]},{μ[3,0],μ[2,1],μ[4,0]-1,μ[2,2]-1,μ[3,1]-μ[1,1]},{μ[1,2],μ[0,3],μ[2,2]-1,μ[0,4]-1,μ[1,3]-μ[1,1]},{μ[2,1],μ[1,2],μ[3,1]-μ[1,1],μ[1,3]-μ[1,1],μ[2,2]-((μ[1,1])^2)}}; (* original variance matrix *)
MatrixForm[Σ]
Dg . Σ . Dg//Simplify(* asymptotic variance of (OverHat[ρ]) *)
Dg . Σ . Dg/.ρ->0 (* asymptotic variance under independence *)

Example: Robust estimation of location

Various ψ-functions of M-estimators that correspond to maximum likelihood estimators in different location models (that is, models where only the location parameter is considered)
For normal distribution ψ(x,t) = (x-t) is unbounded in (x-t) - leads to the sample mean - a non-robust estimator
For Cauchy distribution ψ(x,t) is already bounded, and leads to a robust location estimator
For Laplace distribution (double exponential distribution) ψ(x,t) is the sign of (x-t) - a bounded function that leads the sample median - a very robust estimator
We see a pattern - choosing ψ bounded in (x-t) leads to robust estimators, but we pay by a loss of efficiency for common distributions such as the (uncontaminated) normal distribution.

$Assumptions=Element[θ,Reals];
dist= {NormalDistribution[],StudentTDistribution[1],LaplaceDistribution[],LogisticDistribution[]};
ρ[x_,θ_]=Log[PDF[#,x-θ]]&/@dist;
D[ρ[x,θ],θ]
ψ[x_,θ_]=%;
Plot[Evaluate[ψ[x,3]],{x,3-2,3+2},PlotLegends->{"Normal","Cauchy",  "Laplace", "Logistic"}, PlotLabel->"ψ functions"]

Inverse asymptotic variance (Fisher information)

Clear[x]
-D[ψ[x,θ],{θ}]//Simplify
Table[Integrate[%[[j]]PDF[dist[[j]],x-θ],{x,-∞,∞}],{j,1,Length[dist]}]