Sampling Distribution of ML Estimators: Cauchy Example

Jan Vrbik

We show how to use the Edgeworth series to construct an accurate approximation to the sampling distribution of the maximum likelihood estimator of a parameter of a Cauchy distribution. We then demonstrate the accuracy of this approximation, valid even for relatively small samples.

Introduction

Random variables in statistics have many different distributions; one of them is called Cauchy, and has the following probability density function:

f(x)=(σ/π)·(1/((σ^2)+(((x-μ))^2))),

where x can have any real value. The distribution has two parameters μ and σ, which represent its median (the “location” parameter) and semi-interquartile deviation (the “scale” parameter), respectively. This rather unusual distribution has no mean and infinite standard deviation.

The exact parameter values are usually not known, and need to be estimated by repeating the corresponding random experiment independently n times, and converting the information thus gathered into two respective estimates of μ and σ. The best way of doing this is by maximizing the corresponding likelihood function:

(∏_(i=1)^n)(σ/π)·(1/((σ^2)+((((Subscript[X, i])-μ))^2))),

where the (Subscript[X, i]) are the individual observations, or, equivalently, its logarithm

n ln σ-n ln π-(∑_(i=1)^n)ln[(σ^2)+((((Subscript[X, i])-μ))^2)].

The maximization is usually achieved by setting the first derivative (with respect to each parameter) of the previous expression to zero, thus obtaining the so-called normal equations, namely

(∑_(i=1)^n)(((Subscript[X, i])-(OverHat[μ]))/(((OverHat[σ])^2)+((((Subscript[X, i])-(OverHat[μ])))^2)))=0,

and

(n/(OverHat[σ]))=2(OverHat[σ])(∑_(i=1)^n)(1/(((OverHat[σ])^2)+((((Subscript[X, i])-(OverHat[μ])))^2))),

where the implicit solution (these equations cannot be solved explicitly for μ and σ) defines the maximum likelihood (ML) estimators of μ and σ, denoted (OverHat[μ]) and (OverHat[σ]), respectively. Since these estimators are always implicit functions of the n random observations, they are both random variables, with their own distribution called the sampling distribution.

The normal equations can be rewritten more elegantly as

(OverBar[((X-(OverHat[μ]))/(((OverHat[σ])^2)+(((X-(OverHat[μ])))^2)))])=0

and

(OverBar[(((OverHat[σ])^2)/(((OverHat[σ])^2)+(((X-(OverHat[μ])))^2)))])=(( 1 )/2),

where each bar indicates taking the corresponding sample average—X now symbolically represents all n observations.

Approximating Sampling Distribution of (OverHat[σ])

To simplify our task, we assume that the value of μ is known and equal to zero. From now on, we are thus concerned with estimating only the scale parameter σ, based on the following simplified version of ():

(OverBar[(((OverHat[σ])^2)/(((OverHat[σ])^2)+(X^2)))])=(( 1 )/2).

Furthermore, we can introduce the following so-called pivot

(OverHat[Γ])≡((OverHat[σ])/σ),

which reduces (8) to

(OverBar[(((OverHat[Γ])^( 2))/(((OverHat[Γ])^( 2))+((((X/σ)))^2)))])≡(OverBar[(((OverHat[Γ])^( 2))/(((OverHat[Γ])^( 2))+(U^( 2))))])=(( 1 )/2),

where (Subscript[U, 1]),(Subscript[U, 2]),…,(Subscript[U, n]) is a random independent sample from the Cauchy distribution with μ=0 and the scale parameter equal to 1. It is thus obvious that (OverHat[Γ]) (and its sampling distribution) will be free of σ.

We first expand (OverHat[Γ]) as a power series in ε (a “smallness” parameter—(ε^j) indicates that its coefficient is proportional to (n^(-j/2))), thus:

(OverHat[Γ])=1+ε (Subscript[g, 1])+(ε^2) (Subscript[g, 2])+(ε^3) (Subscript[g, 3])+….

The first term is 1 since the corresponding expansion of (OverHat[σ]) must start with σ—this is a general property of a maximum likelihood estimator. We then similarly expand the expression under the bar in ().

XP[q_,k_]:=Series[q,{ε,0,k}]//Normal//Expand;
lhs=XP[((Γ^2)/((Γ^2)+(U^2)))/.{Γ->1+Sum[(Subscript[g, i])(ε^i),{i,3}]},3]

Placing a bar over each term and then replacing every (OverBar[((U^k)/(((1+(U^k)))^m))]) by 𝔼[((U^k)/(((1+(U^2)))^m))]+ε (OverBar[((U^k)/(((1+(U^2)))^m))-𝔼[((U^k)/(((1+(U^2)))^m))]]), where 𝔼 stands for the corresponding expected value (a plain number), lets us solve for (Subscript[g, 1]), (Subscript[g, 2]), and (Subscript[g, 3]) []. Here Γ1 stands for (OverHat[Γ])-1.

L[q_]:=L[q]=Integrate[(q/(1+(U^2))),{U,-Infinity,Infinity}]/π;
Γ1=Module[{lhs,eq,g},
lhs=XP[((Γ^2)/((Γ^2)+(U^2)))/.{Γ->1+{g[1],g[2],g[3]}.{ε,(ε^2),(ε^3)}},3];eq=XP[lhs/.(U^(k_))(((1+(U^2)))^(m_)):>L[((U^k)/(((1+(U^2)))^(-m)))]+ε (OverBar[((U^k)/(((1+(U^2)))^(-m)))-L[((U^k)/(((1+(U^2)))^(-m)))]])
/.(((1+(U^2)))^(-1)):>L[(1/(1+(U^2)))]+ε(OverBar[(1/(1+(U^2)))-L[(1/(1+(U^2)))]]),3];
Do[g[i]=Expand[g[i]/.Solve[Coefficient[eq,(ε^i)]==0,g[i]][[1]]],{i,3}];{g[1],g[2],g[3]}.{ε,(ε^2),(ε^3)}
]

Cumulants

The next step is to compute the expected value, variance, and the next two cumulants of (OverHat[Γ]). We recall the definition of cumulants as coefficients of (t^i)/i! in the following expansion.

Series[Log[1+Sum[(Subscript[μ, i])(t^i)/i!,{i,2,4}]],{t,0,4}]

This tells us that the second and third cumulants are equal to (Subscript[μ, 2]) and (Subscript[μ, 3]), respectively, but the fourth cumulant is (Subscript[μ, 4])-3((Subscript[μ, 2])^2), where (Subscript[μ, i]) denotes the corresponding central moment, that is,

(Subscript[μ, i])≡𝔼[((((OverHat[Γ])-(Subscript[κ, 1])))^i)],

and (Subscript[κ, 1]) is the expected value of (OverHat[Γ]) (which also happens to be the first cumulant).

To find these, we need a procedure for computing the expected value of expressions involving products and powers of all possible sample averages (our bars) in our previous expansion of Γ1. This is trivial for a single bar (its expected value is always zero), and relatively easy for a product of two or three bars (including powers, as a special case), since one can show that

𝔼[(OverBar[(Subscript[f, 1])(U)])·(OverBar[(Subscript[f, 2])(U)])]=((𝔼[(Subscript[f, 1])(U)·(Subscript[f, 2])(U)])/n),

and

𝔼[(OverBar[(Subscript[f, 1])(U)])·(OverBar[(Subscript[f, 2])(U)])·(OverBar[(Subscript[f, 3])(U)])]=((𝔼[(Subscript[f, 1])(U)·(Subscript[f, 2])(U)·(Subscript[f, 3])(U)])/(n^2)).

For example,

𝔼[(OverBar[(-((1 )/8)+((U^( 2))/(((1+(U^( 2))))^2)))])· ((OverBar[(-(( 1 )/2)+(1/(1+(U^( 2)))))])^2)]=(1/(n π))(∫_(-∞)^∞)(1/(1+(U^( 2))))(-(( 1 )/8)+((U^( 2))/(((1+(U^( 2))))^2)))(((-(( 1 )/2)+(1/(1+(U^( 2))))))^2)dU=-(1/(128n)),

and so on.

Things get more complicated when the overall power is four or higher (we need to go up to six). In such a case, we first need to find all partitions of (Subscript[f, 1])(U),(Subscript[f, 2])(U),(Subscript[f, 3])(U),…, where the individual subsets must have at least two elements each. Thus, for example, having four terms, we either take all of them as a single subset or partition them into two pairs, such as {(Subscript[f, 1]),(Subscript[f, 2])},{(Subscript[f, 3]),(Subscript[f, 4])}—there are three ways of doing this. Similarly, having six (Subscript[f, i]), we can either take all of them or split them into two triplets (10 ways), a pair and a quadruplet (15 ways), or three pairs (15 ways).

The contribution of each of these partitions is then computed using the following scheme: we multiply the expected values of the product of the (Subscript[f, i]) in each subset, multiply the answer by n(n-1)(n-2)…(n+1-k), where k is the number of these subsets, and divide by (n^k). Thus, for example the {(Subscript[f, 1]),(Subscript[f, 2])},{(Subscript[f, 3]),(Subscript[f, 4])},{(Subscript[f, 5]),(Subscript[f, 6])}} partition contributes

((n(n-1)(n-2)·𝔼[(Subscript[f, 1])(Subscript[f, 2])]·𝔼[(Subscript[f, 3])(Subscript[f, 4])]·𝔼[(Subscript[f, 5])(Subscript[f, 6])])/(n^3)).

To find, for example, 𝔼[(OverBar[(Subscript[f, 1])])·(OverBar[(Subscript[f, 2])])·(OverBar[(Subscript[f, 3])])·(OverBar[(Subscript[f, 4])])·(OverBar[(Subscript[f, 5])])·(OverBar[(Subscript[f, 6])])], all such contributions from all corresponding partitions must be added together. This is done by the following program.

EV[x_,io_]:=Quiet@CL[Distribute[H[Expand[x]/.(ε^(j_/;j>2io))->0]]/.H->Q/.ε->1,io];
CL[x_,k_]:=Collect[x,n]/.(n^(i_/;i<-k))->0;
H[q_ w_]:=q H[w]/;FreeQ[q,U];
H[q___,w_ t_,s___]:=Apply[H,Sort[{q,w,t,s}]];
H[w___,((q_)^(i_)),t___]:=Apply[H,Sort[{w,Table[q,{i}],t}//Flatten]];
Quiet[Needs["Combinatorica`"]];
Q[q__]:=Module[{k=Length[{q}],r,s}, r=Select[SetPartitions[{q}],
 FreeQ[Map[Length,#],1]&]; s=Apply[Fa,r,{2}];
 r=Apply[Times,s,{1}]Map[P,Map[Length,s]]/(n^k); Total[CL[r,io]]/.Fa->F];
P[i_]:=Product[n+1-j,{j,i}];
F[q__]:=Module[{w={q}/.(OverBar[a_])->a},L[Apply[Times,w]]];

The second argument of EV specifies the accuracy to be achieved (up to (( 1 )/n) when the argument is equal to 1, up to (1/(n^2)) when it is equal to 2, etc.). At this point, ε becomes a redundant parameter and is eliminated by setting it to 1.

We can thus get the expected value, variance, and the next two cumulants of (OverHat[Γ]).

(Subscript[κ, 1])=1+EV[Γ1,1]
(Subscript[κ, 2])=EV[((Γ1)^2),2]-((EV[Γ1,1])^2)
(Subscript[κ, 3])=EV[((Γ1)^3),2]-3EV[((Γ1)^2),1]EV[Γ1,1]
(Subscript[κ, 4])=CL[EV[((Γ1)^4),3]-4EV[((Γ1)^3),2]EV[Γ1,1]+6EV[((Γ1)^2),1]((EV[Γ1,1])^2)-3((Subscript[κ, 2])^2),3]

Edgeworth Expansion

One can show [] that a good approximation to the probability density function of

Z≡(((OverHat[Γ])-(Subscript[κ, 1]))/(Sqrt[(Subscript[κ, 2])]))

is provided by

((exp(-((z^2)/2)))/(Sqrt[2π]))(1+((Subscript[γ, 3])/6)(Subscript[H, 3])(z)+((Subscript[γ, 4])/24)(Subscript[H, 4])(z)+(((Subscript[γ, 3])^2)/72)(Subscript[H, 6])(z)+…),

where (Subscript[γ, 3])≡((Subscript[κ, 3])/((Subscript[κ, 2])^(3/2))) and (Subscript[γ, 4])≡((Subscript[κ, 4])/((Subscript[κ, 2])^2)) are the normalized cumulants, and (Subscript[H, i])(z) are simple polynomials defined by

(Subscript[H, i])(z)≡(((-1))^i)exp(((z^2)/2))·(((d^i)exp(-((z^2)/2)))/(d(z^i))).

These are closely related to Hermite polynomials; the exact correspondence is clear from the first command.

(Subscript[H, i_])[z_]:=HermiteH[i,(z/(Sqrt[2]))](2^(-i/2))//Expand;
pst=((Exp[-(z^2)/2])/(Sqrt[2π]))(1+((Subscript[κ, 3])/(6((Subscript[κ, 2])^(3/2))))(Subscript[H, 3])[z]+((Subscript[κ, 4])/(24((Subscript[κ, 2])^2)))(Subscript[H, 4])[z]+(((Subscript[κ, 3])^2)/(72((Subscript[κ, 2])^3)))(Subscript[H, 6])[z]);

The resulting expression can be easily transformed into a probability density function of (OverHat[Γ]), which can then be plotted for n=5, for example.

pdf=(pst/.z->((s-(Subscript[κ, 1]))/(Sqrt[(Subscript[κ, 2])])))/(Sqrt[(Subscript[κ, 2])]);
Plot[pdf/.n->5//Simplify,{s,0,4.5}]

Comparison with Empirical Distribution

The best way of verifying the accuracy of this approximation is to generate a large collection of random independent samples of size five from a Cauchy distribution (with the location and scale parameters equal to 0 and 1, respectively), for each such sample compute the value of (OverHat[Γ]) by solving (), and then display these (OverHat[Γ])-values in a histogram.

res=Table[0,{10000}];
Do[u=RandomReal[CauchyDistribution[0,1],5];
res[[i]]=Abs[Γ]/.FindRoot[Total[(Γ^2)/((Γ^2)+(u^2))]==(( 5 )/2),{Γ,1}],{i,10000}]
hist=Histogram[res];

We are now in a position to check how well our approximate probability density function, correspondingly rescaled, matches this histogram (which, for large total frequency, represents—up to small statistical fluctuations—the exact distribution of (OverHat[ Γ])).

Show[hist,Plot[2000pdf/.n->5,{s,0,4.5}],PlotRange->{{0,4.5},{0,1600}}]

Even though this is an improvement over the usual normal approximation, it is still far from a reasonably good match (but one should remember that five is an extremely small sample size—things do improve, rather dramatically, as n increases).

Nevertheless, there is one particular technique one can employ to achieve a substantially better agreement, even with the current value of n=5.

Transforming (OverHat[σ])

We can use the same approach and program to find an approximation for the distribution of T((OverHat[σ])), where T is an arbitrary (preferably increasing) function of a non-negative argument. This can be achieved by defining Tσ and computing its first four cumulants.

Tσ=Series[T[(1+Γ1)σ]-T[σ],{ε,0,3}]//Normal

(Subscript[κ, 1])=T[σ]+EV[Tσ,1]
(Subscript[κ, 2])=Collect[EV[((Tσ)^2),2]-((EV[Tσ,1])^2),n,Simplify]
(Subscript[κ, 3])=EV[((Tσ)^3),2]-3EV[((Tσ)^2),1]EV[Tσ,1]//Simplify
(Subscript[κ, 4])=CL[EV[((Tσ)^4),3]-4EV[((Tσ)^3),2]EV[Tσ,1]+6EV[((Tσ)^2),1]((EV[Tσ,1])^2)-3((Subscript[κ, 2])^2), 3]//Simplify

The third cumulant can be thus made equal to zero for any σ by solving the following differential equation.

DSolve[T'[σ]+σ T''[σ]==0,T[σ],σ]

The simplest nontrivial solution is clear.

T[σ_]:=Log[σ]

Now, we have to recompute the first four cumulants of Tσ.

(Subscript[κ, 1])=T[σ]+EV[Tσ,1]
(Subscript[κ, 2])=Collect[EV[((Tσ)^2),2]-((EV[Tσ,1])^2),n,Simplify]
(Subscript[κ, 3])=EV[((Tσ)^3),2]-3EV[((Tσ)^2),1]EV[Tσ,1]//Simplify
(Subscript[κ, 4])=CL[EV[((Tσ)^4),3]-4EV[((Tσ)^3),2]EV[Tσ,1]+6EV[((Tσ)^2),1]((EV[Tσ,1])^2)-3((Subscript[κ, 2])^2),3]//Simplify

Having reduced (Subscript[κ, 3]) to zero not only greatly simplifies the approximate probability density function of

((ln((OverHat[σ]))-ln(σ))/(Sqrt[ ((2 )/n)+(1/(n^2))])),

but it also makes the corresponding approximation substantially more accurate.

This is now the new approximate probability density function.

pst=((Exp[-(z^2)/2])/(Sqrt[2π]))Series[1+((Subscript[κ, 4])/(24((Subscript[κ, 2])^2)))(Subscript[H, 4])[z],{n,Infinity,1}]//Normal

The improvement can be seen by converting the previous probability density function of () to that of (OverHat[Γ]).

pdf=(pst/.z->((Log[s])/(Sqrt[(Subscript[κ, 2])])))/(Sqrt[(Subscript[κ, 2])])/s;

Let us now compare the last function to the original histogram.

Show[hist,Plot[2000 pdf/.n->5,{s,0,4.5}],PlotRange->{{0,4.5},{0,1600}}]

This has now produced a rather respectable agreement between the new approximation and the exact distribution.

Potential Extensions

In a similar manner, one can construct an approximate joint probability density function for both (OverHat[μ]) and (OverHat[σ]). This requires the corresponding generalization of cumulants, Hermite polynomials, and the Edgeworth expansion. Readers with enough statistical expertise should have no difficulty pursuing that direction on their own since, surprisingly, the required modifications of the program presented in this article would not be that substantial. The most difficult part, namely the routine for computing the expected value of products and powers of sample averages, needs no change at all.

The same technique can also be used to find a good approximation to a probability density function of almost any parameter estimator of a specific distribution.

References

D. F. Andrews and J. E. Stafford, Symbolic Computation for Statistical Inference, New York: Oxford University Press, 2000.

F. Y. Edgeworth, “The Law of Error,” Transactions of the Cambridge Philosophical Society, 20, 1905 pp. 36–66 and 113–141. books.google.com/books/about/Transactions_of_the _Cambridge _Philosophi.html?id=uTs8AAAAMAAJ.