NMST 434, Exercise session VII: M-estimators II
April 2, 2020

Example: M-estimators of location for univariate distributions

Truedist is the true distribution from the location family whose parameter we estimate
Mdist is the distribution whose maximum likelihood estimating equations we use for the construction of M-estimators. In this case, function ψ0 is taken to be -∂Log[f[x;μ]]/∂ μ.
alternatively, function ψ0 can be specified directly as in the case of the Huber loss function or the Tukey biweight function.
In the special case Mdist=Truedist we get the usual maximum likelihood estimators as the M-estimators, and the asymptotic variance in the final computation is, in fact, the Rao-Cramér bound.

Clear["Global`*"]
$Assumptions=Element[μ,Reals]&&b>0;
(*Mdist=NormalDistribution[];*)
Mdist=StudentTDistribution[1];
Truedist =NormalDistribution[];
(*Truedist=LogisticDistribution[];*)
(* Truedist=ExponentialDistribution[1]; *)
(*Truedist=StudentTDistribution[6];*)
f[x_]=PDF[Mdist,x]; (* probability density function of distribution Mdist *)
ψ0[t_]=-((D[f[t],t])/(f[t])); (* ψ function that corresponds to MLE for Mdist *)
(* b=1; *)
ψ0[t_]=Min[b,Max[t,-b]];  (*ψ function for Huber's estimator of location with parameter b *)
(* !!! For non-differentiable functions ψ0, the asymptotic variance matrices computed below may be incorrect - formally, one would have to use Theorem 10 !!! *) 
(*ψ0[t_]=t(((1-((t^2)/(b^2))))^2)Boole[Abs[t]≤b];*) (* ψ function for Tukey's biweight estimator of location with parameter b *)
(*Plot[ψ0[t],{t,-3b/2,3b/2},PlotLabel->"Function ψ"]*)
ψ[x_]=ψ0[x-μ];
Γ=Integrate[(D[ψ[x],{μ}]/.μ->0)PDF[Truedist,x],{x,-∞,∞}] (* Γ matrix: (1,1)-matrix, that is a scalar, now *)
Σ=Integrate[(((ψ[x])^2)/.μ->0)PDF[Truedist,x],{x,-∞,∞}](* Σ matrix: (1,1)-matrix, that is a scalar, now *)
Var1 = Σ/((Γ^2)) (* asymptotic variance of the M-estimator *)

For simple distributions (normal or exponential) we plot the mean (red line), median (orange line), and the estimating function of the M-estimator (as a function of the tuning parameter b). The root of the estimating function is the M-estimate.

OF=Integrate[(ψ[x])PDF[Truedist,x],{x,-∞,∞}]; (* The expectation of the objective function with respect to the true distribution of X, i.e. E ψ[X-μ]. The root of this function is the location parameter that is estimated by the corresponding M-estimator. This location estimator may be different from the mean, or the median, if the distribution is not symmetric *)
(* change Truedist to exponential distribution for an example where the three location estimators are different *)
Manipulate[Plot[OF/.b->B,{μ,-1,2},PlotLabel->"Identifying function",Epilog-> {Directive[Red,Thick],Line[{{Mean[Truedist],-5},{Mean[Truedist],5}}],Directive[Orange,Thick],Line[{{Median[Truedist],-5},{Median[Truedist],5}}]}],{{B,1},0,5}]
Mean[Truedist]
Median[Truedist]
(* to see how the location parameter changes with B = b choose for instance the Huber's ψ function above and re-run the script *)

Computation of the asymptotic variance of Huber’s M-estimator of location using Theorem 10, see Example 40.

Var2 =((Integrate[(x^2)PDF[Truedist,x],{x,-b,b}]+(b^2)(1-CDF[Truedist,b]+CDF[Truedist,-b]))/(((CDF[Truedist,b]-CDF[Truedist,-b]))^2))
%//N

Comparison of the asymptotic variances for Huber’s estimator computed using Theorems 9 and 10. Neither variance expression is true if the distribution is not symmetric, i.e. when μ is misspecified. Run this with Huber’s ψ function chosen above, without specified b parameter.

Plot[{Var1,Var2},{b,0,5},PlotRange->All,PlotStyle->{{Thick},Dashed},PlotLabel->"AVar of Huber's estimator as a function of the parameter b"]

Example: M-estimation of location and scale for univariate distributions

Truedist is the true distribution from the location-scale family whose parameters we estimate
Mdist is the distribution whose maximum likelihood estimating equations we use for the construction of M estimators. In this case, function ψ0 is taken to be -∂Log[f[x;μ,σ]]/∂ μ.
alternatively, function ψ0 can be specified directly as in the case of the Huber loss function or the Tukey biweight function.
In the special case Mdist=Truedist we get the usual maximum likelihood estimators as the M-estimators, and the asymptotic variance in the final computation is, in fact, the Rao-Cramér bound.

Clear["Global`*"]
$Assumptions=Element[μ,Reals]&&σ>0;
(*Mdist=NormalDistribution[];*)
Mdist=StudentTDistribution[1];
Truedist =NormalDistribution[];(*LogisticDistribution[];*)
(*Truedist=StudentTDistribution[6];*)
f[x_]=PDF[Mdist,x];
ψ0[t_]=-((D[f[t],t])/(f[t]));
b=1;
(*ψ0[t_]=Min[b,Max[t,-b]];(* !!! For non-differentiable functions ψ0, the asymptotic variance matrices may be incorrect - formally, one would have to use Theorem 10 !!! *) *)
(*ψ0[t_]=t(((1-((t^2)/(b^2))))^2)Boole[Abs[t]≤b];*)
Plot[ψ0[t],{t,-3b/2,3b/2},PlotLabel->"Function ψ"]
κ=Integrate[ψ0[x]x PDF[Truedist,x],{x,-∞,∞}]; (* standardization constant that ensures correct identification of parameter σ *)
N[κ]
ψ[x_]={ψ0[((x-μ)/σ)],ψ0[((x-μ)/σ)]((x-μ)/σ)-κ};
D[ψ[x],{{μ,σ}}]
Outer[Times,ψ[x],ψ[x]]

Check that the estimator of scale is (Fisher) consistent, i.e. that it identifies σ

ϵ=1/1000;
Integrate[(ψ0[((x-μ)/σ)]((x-μ)/σ)-κ/.μ->0//Simplify)PDF[Truedist,x],{x,-∞,∞}]
Plot[%,{σ,ϵ,Max[b,2]}]
FindRoot[%%,{σ,ϵ,2}]

Matrices Γ and Σ for the computation of influence functions and asymptotic variances

Γ=Integrate[(D[ψ[x],{{μ,σ}}]/.{μ->0,σ->1})PDF[Truedist,x],{x,-∞,∞}]
Σ=Integrate[(Outer[Times,ψ[x],ψ[x]]/.{μ->0,σ->1})PDF[Truedist,x],{x,-∞,∞}]

Influence functions of the M-estimators of location and scale

Plot[Evaluate[-Inverse[Γ].ψ[x]/.{μ->0,σ->1}//Simplify],{x,-3b/2,3b/2},PlotLabel->"Influence functions of the M-estimators of location and scale",PlotLegends->{"Location","Scale"}]
Limit[ψ[t],t->∞]

Asymptotic variance matrix of the M-estimators of location and scale

Inverse[Γ].Σ.Inverse[Γ]/.{μ->0,σ->1}//Simplify//MatrixForm
%//N//MatrixForm

Fisher information matrix of the true distribution, to compare with the asymptotic variance of the M-estimators above

U=-((D[(1/σ)PDF[Truedist,((t-μ)/σ)],{{μ,σ}}])/((1/σ)PDF[Truedist,((t-μ)/σ)]));
FIM = Integrate[Outer[Times,U,U]PDF[Truedist,t]/.{μ->0,σ->1},{t,-∞,∞}];
Inverse[FIM]//MatrixForm
%//N//MatrixForm

Derivative of log(det(Σ)) is ((((Σ^(-1))))^T)

d=5; (* dimension *)
A=D[Log[Det[Table[(Subscript[x, i,j]),{i,1,d},{j,1,d}]]],{Table[(Subscript[x, i,j]),{i,1,d},{j,1,d}]}];
Short[A]
B=Transpose[Inverse[Table[(Subscript[x, i,j]),{i,1,d},{j,1,d}]]];
Cond=Flatten[Table[(Subscript[x, i,j])->RandomReal[],{i,1,d},{j,1,d}]] (* Evaluate the matrices randomly *)
(A/.Cond)
(B/.Cond)

Example: Multivariate M-estimation of location and scatter

Clear["Global`*"];
$Assumptions=μ1 ∈Reals && μ2 ∈Reals&&σ1>0 && σ2>0 && -1<ρ<1;
μ={μ1,μ2};
Σ={{((σ1)^2),ρ σ1 σ2},{ρ σ1 σ2, ((σ2)^2)}};
θ={μ1,μ2,σ1,ρ,σ2};
Truedist =MultinormalDistribution[DiagonalMatrix[{1,1}]];
b=1;
ψ0[t_]=Min[b,Max[t,-b]]/t;(* !!! For non-differentiable functions ψ0, the asymptotic variance matrices may be incorrect - formally, one would have to use Theorem 10 !!! *) 
(*ψ0[t_]=1;*)
(*ψ0[t_]=(((2+1))/((t^2)+1));*)(* Corresponds to t-distribution with 1 degree of freedom - Cauchy distribution *)
(*ψ0[t_]=t(((1-((t^2)/(b^2))))^2)Boole[Abs[t]≤b]/t;*)
Plot3D[Evaluate[ψ0[√(({x1,x2}-{0,0}).Inverse[DiagonalMatrix[{1,1}]].({x1,x2}-{0,0}))]{x1,x2}],{x1,-3b/2,3b/2},{x2,-3b/2,3b/2},PlotRange->All,PlotLabel->"Functions ψ for location parameters",ImageSize->Large]
Plot3D[Evaluate[ψ0[√(({x1,x2}-{0,0}).Inverse[DiagonalMatrix[{1,1}]].({x1,x2}-{0,0}))]{(x1^2),x1 x2,(x2^2)}],{x1,-3b/2,3b/2},{x2,-3b/2,3b/2},PlotRange->All,PlotLabel->"Functions ψ for scatter parameters",ImageSize->Large]
κ=NIntegrate[ψ0[√(({x1,x2}-{0,0}).Inverse[DiagonalMatrix[{1,1}]].({x1,x2}-{0,0}))](x1^2)PDF[Truedist,{x1,x2}],{x1,-∞,∞},{x2,-∞,∞}]
ψ[x1_,x2_]=Delete[Join[ψ0[√(({x1,x2}-μ).Inverse[Σ].({x1,x2}-μ))]({x1,x2}-μ),Flatten[((ψ0[√(({x1,x2}-μ).Inverse[Σ].({x1,x2}-μ))])/κ)Outer[Times,{x1,x2}-μ,{x1,x2}-μ]-Σ]],5];

Check that we correctly identified all the parameters. Integrals should be all zero at the true value of the parameters.
(addition and subtraction of 1 only to avoid warnings “Integral may be zero”)

NIntegrate[(1+(ψ[x1,x2]/.{μ1->0,μ2->0,σ1->1,σ2->1,ρ->0}))PDF[Truedist,{x1,x2}],{x1,-∞,∞},{x2,-∞,∞},AccuracyGoal->10,WorkingPrecision->10]-1

If the integrals in Γ and V can be computed exactly (only for the estimating equations from the multivariate normal distribution)

Γ=Integrate[(D[ψ[x1,x2],{{μ1,μ2,σ1,ρ,σ2}}]/.{μ1->0,μ2->0,σ1->1,σ2->1,ρ->0})PDF[Truedist,{x1,x2}],{x1,-∞,∞},{x2,-∞,∞}]
V=Integrate[(Outer[Times,ψ[x1,x2],ψ[x1,x2]]/.{μ1->0,μ2->0,σ1->1,σ2->1,ρ->0})PDF[Truedist,{x1,x2}],{x1,-∞,∞},{x2,-∞,∞}]

Numerical computation of matrices Γ and V

Γ=NIntegrate[(1+D[ψ[x1,x2],{{μ1,μ2,σ1,ρ,σ2}}]/.{μ1->0,μ2->0,σ1->1,σ2->1,ρ->0})PDF[Truedist,{x1,x2}],{x1,-∞,∞},{x2,-∞,∞},AccuracyGoal->10,WorkingPrecision->10]-1
V=NIntegrate[(1+Outer[Times,ψ[x1,x2],ψ[x1,x2]]/.{μ1->0,μ2->0,σ1->1,σ2->1,ρ->0})PDF[Truedist,{x1,x2}],{x1,-∞,∞},{x2,-∞,∞},AccuracyGoal->10]-1

Influence function of the M-estimators of location and scatter

IF=Evaluate[-Inverse[Γ].ψ[x1,x2]/.{μ1->0,μ2->0,σ1->1,σ2->1,ρ->0}//Simplify];
Plot3D[IF,{x1,-3b/2,3b/2},{x2,-3b/2,3b/2},PlotLabel->"Influence functions of bivariate M-estimators of location and scatter",PlotLegends->{"(Subscript[μ, 1])","(Subscript[μ, 2])","(Subscript[σ, 1])","ρ","(Subscript[σ, 2])"},AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])","IF((Subscript[x, 1]),(Subscript[x, 2]))"},ImageSize->Large,PlotRange->All]
Plot3D[IF[[4]],{x1,-3b/2,3b/2},{x2,-3b/2,3b/2},PlotLabel->"Influence function of the M-estimator of correlation coef.",PlotLegends->{"ρ"},AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])","IF((Subscript[x, 1]),(Subscript[x, 2]))"},ImageSize->Large,PlotRange->All]

Asymptotic variance matrix of the M-estimator of location and scatter

Chop[Inverse[Γ].V.Inverse[Γ],(10^(-10))]//MatrixForm

Inverse Fisher information matrix for multivariate normal distribution

U=-((D[PDF[MultinormalDistribution[{μ1,μ2},Σ],{x1,x2}],{θ}])/(PDF[MultinormalDistribution[{μ1,μ2},Σ],{x1,x2}]))//Simplify;
FIM = Integrate[(Outer[Times,U,U]PDF[MultinormalDistribution[DiagonalMatrix[{1,1}]],{x1,x2}])/.{μ1->0,μ2->0,σ1->1,σ2->1,ρ->0},{x1,-∞,∞},{x2,-∞,∞}];
Inverse[FIM]//MatrixForm