Asymptotic Expansions for Moments of Functions of Random Sequences

http://library.wolfram.com/infocenter/Conferences/7454/

Jan Hurt

Charles University of Prague, Czech Republic

Introduction

We consider a smooth real function g and a random sequence (Subscript[{(Subscript[T, n])},  n=1,...,∞]) which converges to a parameter θ in some sense. In Statistics, this sequence is often a sequence of statistics that estimate some unknown parameter based on random sample of size n. We give the asymptotic expansions (theoretically of arbitrary order) of the expected value, variance, and mean squared error of the sequence (Subscript[{g((Subscript[T, n]))} , n=1,...,∞]). The method is also known as the method of statistical differentials or delta method. A supporting theory may be found in [1], [2], and [3] e. g.

Asymptotics in one-dimensional case

The expectation

If (Subscript[{g(t,n)}, n=1,...∞]) is a sequence of functions that for some q>0 admits a continuous (((q+1))^st) derivative with respect to the first argument in some neighbourhood of θ and E|(Subscript[T, n])-θ((|)^(q+1))⟶0 as n⟶∞ then under general conditions

E g((Subscript[T, n]),n)=g(θ,n)+(∑_(j=1)^q)(1/(j!))(g^((j)))(θ,n)((E((Subscript[T, n])-θ))^j)+O(E|(Subscript[T, n])-θ((|)^(q+1))).

The derivatives of g are taken with respect to the first argument.

The variance

If E|(Subscript[T, n])-θ((|)^(q+2))⟶0 as n⟶∞ then under general conditions

var g((Subscript[T, n]),n)=(Underscript[(∑_(j=1)^q)(∑_(k=1)^q), j+k≤ q+1])(1/(j!))(1/(k!))(g^((j)))(θ,n)(g^((k)))(θ,n)cov[((((Subscript[T, n])-θ))^j),((((Subscript[T, n])-θ))^k)]+O(E|(Subscript[T, n])-θ((|)^(q+2))).

Asymptotics in multidimensional case

The expectation

Let (Subscript[{(Subscript[T, n])},  n=1,...,∞])=(Subscript[{(Subscript[T, n ,1]),...(Subscript[T, n,r])}, n=1,...,∞]) be an r-dimensional sequence of random vectors and let (Subscript[{g(t,n)}, n=1,...∞]) be a sequence of real functions defined on (Subscript[R, r])×𝒩 which is (q+1) times totally differentiable with respect to t in some neighbourhood of θ and (max_(1≤i≤r))E|(Subscript[T, n,i])-(Subscript[θ, i])((|)^(q+1))⟶0 as n⟶∞ then under general conditions

E g((Subscript[T, n]),n)=g(θ,n)+(∑_(j=1)^q)(1/(j!))(Underscript[∑...∑, (Subscript[m, 1])+...+(Subscript[m, r])=j])(D^j)g((Subscript[m, 1]),...,(Subscript[m, r]);θ)E[((((Subscript[T, n,1])-(Subscript[θ, 1])))^(Subscript[m, 1]))...((((Subscript[T, n,r])-(Subscript[θ, r])))^(Subscript[m, r]))]+O((max_(1≤i≤r))E|(Subscript[T, n,i])-(Subscript[θ, i])((|)^(q+1)))

where (Subscript[m, 1]),...,(Subscript[m, r]) are all compositions of integer j into r parts, and (D^j)g((Subscript[m, 1]),...,(Subscript[m, r]);θ)=(Subscript[((∂^j)g/∂(Subscript[y, 1]^(Subscript[m, 1]))...∂(Subscript[y, r]^(Subscript[m, r]))), y=θ]) .   We will see that the formula is more transparent in Mathematica language.

The variance

Similarly, if (max_(1≤i≤r))E|(Subscript[T, n,i])-(Subscript[θ, i])((|)^(q+2))⟶0 as n⟶∞ then under general conditions

var g((Subscript[T, n]),n)=(Underscript[(∑_(j=1)^q)(∑_(k=1)^q), j+k≤q+1  ])(1/(j!))(1/(k!))(Underscript[∑...∑, (Subscript[m, 1])+...+(Subscript[m, r])=j])   (Underscript[∑...∑, (Subscript[n, 1])+...+(Subscript[n, r])=k]) (D^j)g((Subscript[m, 1]),...,(Subscript[m, r]);θ)(D^k)g((Subscript[n, 1]),...,(Subscript[n, r]);θ)×cov[((((Subscript[T, n,1])-(Subscript[θ, 1])))^(Subscript[m, 1]))...((((Subscript[T, n,r])-(Subscript[θ, r])))^(Subscript[m, r])),((((Subscript[T, n,1])-(Subscript[θ, 1])))^(Subscript[n, 1]))...((((Subscript[T, n,r])-(Subscript[θ, r])))^(Subscript[n, r]))]+O((max_(1≤i≤r))E|(Subscript[T, n,i])-(Subscript[θ, i])((|)^(q+2)))

where (Subscript[n, 1]),...,(Subscript[n, r]) are all compositions of integer j into r parts.

Mean squared error

Sometimes, an asymptotic expansion of the mean squared error (MSE) is needed rather than that of variance. If w(θ) is a parametric function in question then a recommendable formula for calculating MSE is

((E [g((Subscript[T, n]),n)-w(θ)])^2)=var g((Subscript[T, n]),θ)+(([Eg,((Subscript[T, n]),θ)-w(θ)])^2).

Sequences of statistics

If (Subscript[{(Subscript[T, n])}, n=1,...,∞]) are statistics based on a random sample of size n, then we often have

E|(Subscript[T, n])-θ((|)^s)=((O(n)^(-δ,s)))

for some δ>0. For example, δ=1/2 if (Subscript[T, n]) is a 2(k^th) sample central moment from a parent population with finite 2k,(s^th) moment or δ=1 for some extreme statistics.

One-dimensional case

Expected value

In the Mathematica procedure, function g should be defined as a function of two arguments, g:=g(t,n) even in case that g does not depend on n explicitly, the power density function of (Subscript[T, n]) as a function of three arguments, p,d,f:=p,d,f(x,θ,n). Parameter range is the support of the corresponding density in the form of a list. (This parameter may be considered as superfluous at the first look. Unfortunately, sometimes for some densities, even correctly defined on the whole interval (-∞,∞), the integration part of the procedure does not work.). Parameter assumptions is also used in the integration part.

Implementation

Clear[expected];
expected[g_,θ_,n_,q_,pdf_,range_,assumptions___]:=
Module[{j,t,x,derivatives,centralmoments,ass,aux},
ass=assumptions&&n>0&&q>0//Simplify;
derivatives=Table[D[g[t,n],{t,j}]/.t->θ,{j,2q}]//Simplify;
centralmoments=Table[Integrate[Expand[(((x-θ))^j)]pdf[x,θ,n],{x,range}//Flatten,Assumptions->ass]//Simplify,{j,2q}];
aux=g[θ,n]+(1/(Range[2q]!)).(derivatives*centralmoments)//Simplify;
Series[aux,{n,∞,q}//Simplify]
]

Example 1 (Estimation of reliability in the exponential case)

We consider an exponential parent distribution Exp(θ), with the distribution function F(x;θ)=1-(e^(-x/θ)), the reliability (or survival function) R(x;θ)=1-F(x;θ)=(e^(-x/θ)), and the density f(x;θ)=(1/θ)(e^(-x/θ)) for x>0 and zero otherwise. We will estimate R:=R(x;θ) at some fixed mission time x. Let (OverBar[X]) be the mean from the sample of size n. The best unbiased estimate of R (see [4]) is

(Subscript[R, 1])=(((1-(x/(n (OverBar[X])))))^(n-1))

for (OverBar[X])>x/n and 1 otherwise. (This will not be emphasized in the sequel.) The maximum likelihood estimate is

(Subscript[R, 2])=(e^(-x/(OverBar[X]))).

Suppose now that parameter λ=1/θ is a random variable with the a priori distribution of the gamma-type with parameters p and 1/α: GammaDistribution[p,1/α]. Then the a posteriori distribution is also of the gamma type with parameters p+n and 1/(α+∑(Subscript[X, i])) with the density

PDF[GammaDistribution[p+n,1/(α+Σ)],λ]

where Σ stands for (Subscript[ΣX, i]). The Bayesian estimate of the reliability is the expected value of R=(e^(-λ,x)) with respect to the posterior density of λ.

r3=Assuming[n>0&&p>0&&Σ>0&&α>0&&x>0,(∫_0^∞)(E^(-λ x))PDF[GammaDistribution[p+n,1/(α+Σ)],λ]\[DifferentialD]λ]

so that the Bayesian estimate is

(Subscript[R, 3])=(((((α+n (OverBar[X]))/(α+x+n (OverBar[X])))))^(n+p)).

Sample mean (OverBar[X]) possesses the density

(pdfgamma[x_,θ_,n_]=PDF[GammaDistribution[n,θ/n],x])//TraditionalForm

Define functions g needed in the asymptotic expansions for the above three estimates. BUE, MLE, and Bayes are the abbreviations for the best unbiased, maximum likelihood, and Bayesian estimates, respectively.

Clear[x,α,p,gexpBUE,gexpMLE,gexpBayes];
gexpBUE[t_,n_]:=(((1-(x/(t n))))^(n-1))
gexpMLE[t_,n_]:=Exp[-x/t]
gexpBayes[t_,n_]:=(((((α+t n)/(α+x+t n))))^(n+p))

As we will see later, the expansions are the functions of the ratio κ=x/θ only.

Obviously, the asymptotic expansion of the expected value of the unbiased estimate is unbiased up to the desired order:

(expectedBUE=expected[gexpBUE,θ,n,2,pdfgamma,{0,∞},n>0&&x>0&&θ>0])//TraditionalForm

expected[gexpBUE,θ,n,5,pdfgamma,{0,∞},n>0&&x>0&&θ>0]//TraditionalForm

For the maximum likelihood estimate we have

(expectedMLE=expected[gexpMLE,θ,n,1,pdfgamma,{0,∞},n>0&&x>0&&θ>0])//TraditionalForm

and using the transform κ=x/θ we have the normal form

(expectedMLEκ=(E^(-(x/θ)))-((x (E^(-(x/θ))) (2 θ-x))/(2 n (θ^2)))/.{x->θ κ,x/θ->κ}//Cancel)//TraditionalForm

For the Bayesian estimate

(expectedBayes=expected[gexpBayes,θ,n,1,pdfgamma,{0,∞},n>0&&x>0&&θ>0&&α>0])//TraditionalForm

If we put β=α/θ we have

(expectedBayesκ=(E^(-(x/θ)))-((x (E^(-(x/θ))) (p θ-x-α+θ))/(n (θ^2)))/.{x->θ κ,x/θ->κ,α->θ β}//Expand//Cancel//Collect[#,(E^(-κ))]&//Collect[#,n]&)//TraditionalForm

Next we give the graphical comparison of the bias of the maximum likelihood and Bayesian estimates. Note that the shape of the curves does not depend on n.

Manipulate[Plot[{(((E^(-κ)) (-2+κ) κ)/(2 n)),(((E^(-κ)) (-κ-p κ+β κ+(κ^2)))/n)},{κ,0,5},PlotRange->All,AxesLabel->{"κ","bias"}],{{n,20},{5,10,20,50,100,1000}},{{p,5},1,20},{{β,1},0.1,10}]

Example 2 (A simple system reliability)

Suppose that (Subscript[X, 1]),...,(Subscript[X, n]) is a random sample from the uniform distribution on (0,θ), θ an unknown parameter. The (Subscript[X, (n)])=max ((Subscript[X, 1]),...,(Subscript[X, n])) possesses the density function

Clear[pdfmax];
pdfmax[x_,θ_,n_]:=((n (x^(n-1)))/(θ^n))

for 0<x<θ and zero otherwise. For x>0 fixed in the sequel we need to estimate the parametric function w(θ,x)=exp(-2x/θ). This function arises in reliability theory. A system consists of two parts with the same average lifetime θ/2. Suppose that the distribution of the time to failure of the first part is exponential and that of the second part is uniform on (0,θ). If the second part fails down, the failure of the system is instantly observable. If the second part fails down, however, the system operates mistakenly but its wrong function is not immediately observable. We want to estimate the reliability of the first part assuming that the observations of the times to failure (Subscript[X, 1]),...,(Subscript[X, n]) of the second part are available. The best unbiased estimate of θ is (OverHat[θ]) = (n+1)(Subscript[X, m,a,x])/n where (Subscript[X, m,a,x]) is the maximum of the observations so it is natural to estimate w by

(OverHat[w])(θ,x)=exp(-((2x,n)/(n+1))(1/(Subscript[X, m,a,x]))).

Clear[x,w];
w[t_,n_]:=Exp[-((2x n)/(n+1))(1/t)]

expected[w,θ,n,2,pdfmax,{0,θ},n>0&&x>0&&θ>0]//TraditionalForm

Alternatively expressed in terms of κ:

(E^(-((2 x)/θ)))-((2 (x (E^(-((2 x)/θ))) (θ-x)))/((n^2) (θ^2)))/.{x->θ κ,x/θ->κ}//Expand//Cancel//Simplify//Collect[#,n]&//TraditionalForm

We see that the estimate is unbiased up to the second order.

Note that the direct approach does not work:

expectationdirect=Assuming[n∈Integers&&n>0&&θ>0&&x>0,(∫_0^θ)Exp[-((2x n)/(n+1))(1/t)]pdfmax[t,θ,n]\[DifferentialD]t]//FunctionExpand//Simplify
Series[expectationdirect,{n,∞,2}]//Simplify

Variance

Implementation

Clear[variance];
variance[g_,θ_,n_,q_,pdf_,range_,assumptions___]:=
Module[{j,u,v,t,x,dd,centralmoments4q,centralmoments2q,covariances,ass,aux},
ass=assumptions&&n>0&&q>0//Simplify;
dd=Table[D[g[t,n],{t,j}]/.t->θ,{j,2q}]//Simplify;
centralmoments4q=Table[Integrate[Expand[(((x-θ))^j)]pdf[x,θ,n],{x,range}//Flatten,Assumptions->ass]//Simplify,{j,4q}];
centralmoments2q=Take[centralmoments4q,2q];
covariances=Partition[Rest[centralmoments4q],2q,1]-Outer[Times,centralmoments2q,centralmoments2q];
aux=((1/(Range[2q]!))*dd).covariances.((1/(Range[2q]!))*dd);
Series[aux,{n,∞,q}//Simplify]
]

Example 1 (continued)

The best unbiased estimate:

varBUE=variance[gexpBUE,θ,n,2,pdfgamma,{0,∞},n>0&&x>0&&θ>0]

therefore

varBUEκ=varBUE/.{x->θ κ,x/θ->κ}//Simplify//TraditionalForm

The maximum likelihood estimate:

varMLE=variance[gexpMLE,θ,n,2,pdfgamma,{0,∞},n>0&&x>0&&θ>0]

varMLEκ=varMLE/.{x->θ κ,x/θ->κ}//Simplify//TraditionalForm

The Bayesian estimate:

varBayes=variance[gexpBayes,θ,n,2,pdfgamma,{0,∞},n>0&&x>0&&θ>0]

varBayesκ=varBayes/.{x->θ κ,x/θ->κ,α->θ β}//Simplify//TraditionalForm

An alternative approach

If we put h((Subscript[T, n]),n):=(([g((Subscript[T, n]),n)-E,g((Subscript[T, n]),n)])^2) then from the definition of the variance it follows that

var g((Subscript[T, n]),n)=((E[g((Subscript[T, n]),n)-E g((Subscript[T, n]),n)])^2)=E,h((Subscript[T, n]),n)

so that we may carefully use the formula for the asymptotic expansion of the expected value.

Clear[varianceAlternative];
varianceAlternative[g_,θ_,n_,q_,pdf_,range_,assumptions___]:=Module[{h,t,n1},
h[t_,n1_]:=(((g[t,n1]-expected[g,θ,n1,q,pdf,range,assumptions]//Normal))^2)//Simplify;
expected[h,θ,n,q,pdf,range,assumptions]
	]

Some illustrations:

varianceAlternative[gexpBUE,θ,n,2,pdfgamma,{0,∞},n>0&&x>0&&θ>0]

varianceAlternative[gexpMLE,θ,n,2,pdfgamma,{0,∞},n>0&&x>0&&θ>0]

varianceAlternative[gexpBayes,θ,n,2,pdfgamma,{0,∞},n>0&&x>0&&θ>0]

Next we give the comparison of the timings for doing a computation in the two different ways.

{Timing[variance[gexpBUE,θ,n,2,pdfgamma,{0,∞},n>0&&x>0&&θ>0];],Timing[varianceAlternative[gexpBUE,θ,n,2,pdfgamma,{0,∞},n>0&&x>0&&θ>0];]}

Multidimensional case

Expected value

Implementation

Here we effectively use the Compositions[] function from Combinatorica package.

Clear[expectedMulti];
Needs["Combinatorica`"];
expectedMulti[g_,θ_,n_,q_,pdf_,range_,assumptions___]:=
Module[{j,tt,t,x,r,compositions,derivatives,centralmoments,ass,aux,ttrange,integrallimits},
If[Not[Length[θ]==Length[range]],Print["Incompatible dimenensions of the input vectors."];Return[]];
ass=assumptions&&n>1&&q>0//Simplify;
r=Length[θ];
compositions=Table[Compositions[j,r],{j,2q}];
tt=Array[t,r];
derivatives=(Table[Map[D[g[tt,n],Apply[Sequence,Transpose[{tt,#}]]]&,compositions[[j]]],{j,2q}]//Simplify)/.Thread[tt->θ];
ttrange={tt,Transpose[range]};
integrallimits=Transpose[Partition[Flatten[ttrange],r]];
integrallimits=Apply[Sequence,integrallimits];
centralmoments=Table[Map[Integrate[# pdf[tt,θ,n],integrallimits,Assumptions->ass]&,Map[Expand[Apply[Times,#]]&,Map[(((tt-θ))^(#))&,compositions[[j]]]]]//Simplify,{j,2q}];
aux=Flatten[derivatives,1].Flatten[(1/(Range[2q]))*centralmoments,1];
aux=Series[g[θ,n]+aux,{n,∞,q}//Simplify]]

Example 4 (Estimation of the CDF of a normal distribution)

Let X be distributed as N(μ,σ) with mean μ and standard deviation σ, both the parameters unknown. For a fixed real c we will estimate the distribution function

P{X<c}=Φ(((c-μ)/σ)).

where Φ stands for CDF of N(0,1).  The statistics (s^2)=(1/(n-1))(∑_(i=1)^n)((((Subscript[X, i])-(OverBar[X])))^2) is an unbiased estimate of (σ^2). The maximum likelihood estimate of P{X<c} is

Φ((c-(OverBar[X]))/(Sqrt[((n-1)/n)(s^2)])) ,

so that the estimating function is

Clear[gnormaldistributionMLE,c];
gnormaldistributionMLE[{t1_,t2_},n_]:=CDF[NormalDistribution[0,1],(c-t1)/(Sqrt[((n-1)/n)t2])]

The join distribution of (OverBar[X]) and (s^2) has the density (here σ2 stands for (σ^2)):

Clear[pdfnormal];
pdfnormal[{x1_,x2_},{μ_,σ2_},n_]:=PDF[NormalDistribution[μ,(Sqrt[σ2/n])],x1]*((n-1)/(σ2))*PDF[ChiSquareDistribution[n-1],((n-1)/(σ2))x2]

The asymptotic expansion of the expected value of the above estimate is therefore

(expectedMLEnormal=expectedMulti[gnormaldistributionMLE,{μ,σ2},n,1,pdfnormal,{{-∞,∞},{0,∞}},n>1&&μ>0&&σ2>0])

Note that the first term is actually Φ((c-μ)/σ). Put κ=(c-μ)/σ. Then after some algebra we get a more concise expression for the above asymptotic expansion:

expectedMLEnormalκ=(expectedMLEnormal/.σ2->(σ^2)//PowerExpand//ExpandAll)/.{c-μ->σ κ,(c-μ)/σ->κ,σ->(c-μ)/κ}//Simplify

Traditional form of the last expression with φ denoting the density function of N(0,1) reads

Φ(κ)+(1/(4n))κ φ(κ)((κ^2)-3)+O((n^(-2))).

Example 5 (Estimation of the distribution function of a truncated exponential distribution if the point of the truncation is not known)

This is an example of an unsuccessful attempt to use the above procedure. Consider a population with the distribution function F(x;θ,a) defined as

Clear[cdfexptruncated];cdfexptruncated[x_,θ_,a_]:=Which[x>a,1-Exp[-(1/(θ-a))(x-a)],
True,0]

with both parameters θ, and a unknown. Here a plays the role of the point of truncation. We want to estimate this function at a fixed mission time x:=c. The mean of the distribution is θ:

mean=Integrate[x D[cdfexptruncated[x,θ,a],x],{x,0,∞},Assumptions->θ>a>0]

It is not difficult to see that the maximum likelihood estimates of θ and a are the sample mean (OverBar[X]) and the minimum of the observations (Subscript[X, (1)])=min((Subscript[X, 1]),...,(Subscript[X, n])), respectively. It is therefore natural to estimate F(c;θ,a) by the statistic

(OverHat[R ])=1-exp{-((c-(Subscript[X, (1)]))/((OverBar[X])-(Subscript[X, (1)])))}.

if (Subscript[X, (1)])<c and 1 otherwise. The join density of (OverBar[X]) and (Subscript[X, (1)]) is (the arguments x1 and x2 in the same order):

Clear[pdfexptruncated];
pdfexptruncated[{x1_,x2_},{θ_,a_},n_]:=Which[a<x2<x1,(1/((n-2)!))((n^n)/(((θ-a))^n))(((x1-x2))^(n-2))Exp[-n((x1-a)/(θ-a))],
	 True,0]

The estimated function is

Clear[gexptruncated,c];
gexptruncated[{t1_,t2_},n_]:=1-Exp[-((c-t2)/(t1-t2))]

We will illustrate the problem on the calculation of the variance of the sample mean using the join density. We observe the dependence on the order of integration. The following integral is easily done:

Integrate[Expand[(((x1-θ))^2)]pdfexptruncated[{x1,x2},{θ,a},n],{x2,0,∞},{x1,0,∞},Assumptions->n>1&&0<a<θ]//Simplify

But it is difficult to get the result with reverse order of integration in a finite time.

Integrate[Expand[(((x1-θ))^2)]pdfexptruncated[{x1,x2},{θ,a},n],{x1,0,∞},{x2,0,∞},Assumptions->n>1&&0<a<θ]//Simplify

An open problem is how to automate the calculation of the moments in such situations. Of course, the procedure cannot give the answer as well:

expectedMulti[gexptruncated,{θ,a},n,1,pdfexptruncated,{{0,∞},{0,∞}},n>1&&θ1>0&&θ2>0&&θ1>θ2]

We can however incorporate the general idea and proceed step by step:

q=2;
Needs["Combinatorica`"];
compositions=Table[Compositions[j,2],{j,2q}];
compositionsflattened=Flatten[compositions,1];
derivatives=(Table[Map[D[gexptruncated[{x1,x2},n],Apply[Sequence,Transpose[{{x1,x2},#}]]]&,compositions[[j]]],{j,2q}]//Simplify)/.Thread[{x1,x2}->{θ,a}];
derivativesmultipiedflattened=Flatten[(1/(Range[2q]!))derivatives,1];
moments=Map[Integrate[Expand[(((x2-a))^(#[[1]])) ]Integrate[Expand[(((x1-θ))^(#[[2]]))]pdfexptruncated[{x1,x2},{θ,a},n],{x1,0,∞},Assumptions->θ>a>0&&n>1&&a<x2],{x2,a,∞},Assumptions->θ>a>0&&n>1]&//Simplify//PowerExpand,compositionsflattened];
aux=derivativesmultipiedflattened.moments;
aux=Series[gexptruncated[{θ,a},n]+aux,{n,∞,q}]//Simplify

Variance

The procedure for the variance is a technical extension for that of the expectation. But it is also possible to utilize the formula for the expected value as we did in the one-dimensional case.

Implementation of an alternative approach

Clear[varianceMultiAlternative];
varianceMultiAlternative[g_,θ_,n_,q_,pdf_,range_,assumptions___]:=Module[{h,t,n1},
h[t_,n1_]:=(((g[t,n1]-expectedMulti[g,θ,n1,q,pdf,range,assumptions]//Normal))^2)//Simplify;
expectedMulti[h,θ,n,q,pdf,range,assumptions]
	]

Example 4 (continued)

Warning! The following calculation takes some 3/4 of an hour on an average PC!

varMLEnormal=varianceMultiAlternative[gnormaldistributionMLE,{μ,σ2},n,2,pdfnormal,{{-∞,∞},{0,∞}},n>1&&n∈Integers&&σ2>0]

varMLEnormalκ=(varMLEnormal/.σ2->(σ^2)//PowerExpand//ExpandAll)/.{c-μ->σ κ,(c-μ)/σ->κ,σ->(c-μ)/κ}//Simplify

Deficiency

In order to compare two equally efficient statistical procedures Hodges and Lehmann [5] introduced the concept of deficiency. Suppose we have two statistical procedures A and B, say, and their efficiencies are measured by quantities (Subscript[V, n]^A) and (Subscript[V, n]^B), respectively, where n is the size of the underlying random sample. Usually V's are mean squared errors. Using formulas given above, these quantities may often be represented as

(Subscript[V, n]^A)=(c/(n^r))+(a/(n^(r+1)))+o((n^(-r-1))),

(Subscript[V, n]^B)=(c/(n^r))+(b/(n^(r+1)))+o((n^(-r-1))).

Note that c must be the same in both expansions since we suppose both procedures equally efficient. Then the (asymptotic) deficiency of B with respect to A is defined as

(Subscript[d, B,A])=((b-a)/(r,c)).

Roughly speaking, to attain the same value of the mean squared error of B and A we need some (Subscript[d, B,A]) more observation for procedure B.

Example 1 (continued)

Mean squared errors for the three estimates:

mseBUE=(varBUE//ExpandAll)/.{x->θ κ,x/θ->κ,α->θ β}//Simplify
mseMLE=Series[varMLE+(((expectedMLE-(E^(-(x/θ)))))^2),{n,∞,2}]/.{x->θ κ,x/θ->κ,α->θ β}//Simplify
mseBayes=Series[varBayes+(((expectedBayes-(E^(-(x/θ)))))^2),{n,∞,2}]/.{x->θ κ,x/θ->κ,α->θ β}//Simplify

Now the deficiencies are calculated by the following scheme using the built-in function SeriesCoefficient. Therefore the deficiency (Subscript[d, M,L,E,B,U,E]) is:

(SeriesCoefficient[mseMLE,2]-SeriesCoefficient[mseBUE,2])/SeriesCoefficient[mseBUE,1]//Simplify

Similarly we can proceed with (Subscript[d, B,a,y,e,s,B,U,E]):

(SeriesCoefficient[mseBayes,2]-SeriesCoefficient[mseBUE,2])/SeriesCoefficient[mseBUE,1]//Simplify

And (Subscript[d, B,a,y,e,s,M,L,E]) is

(SeriesCoefficient[mseBayes,2]-SeriesCoefficient[mseMLE,2])/SeriesCoefficient[mseBUE,1]//Simplify

Manipulate[Plot[{3-5 κ+((5 (κ^2))/4),3+(p^2)+(β^2)-8 κ+3 (κ^2)-2 p (-2+β+2 κ)+β (-6+4 κ),(p^2)-6 β+(β^2)-3 κ+4 β κ+((7 (κ^2))/4)-2 p (-2+β+2 κ)},{κ,0,5},PlotRange->All,PlotStyle->{Red,Green,Blue},AxesLabel->{"κ","deficiency"}],{{p,5},1,20},{{β,1},0.1,10}]

Acknowledgments

The research has been supported by Grant of Czech Ministry of Education MSM 0021620839, Methods of modern mathematics and their applications.

References

[1] Cramer, H. (1946): Mathematical Methods of Statistics. Princeton, Princeton University Press. 
[2] Prakasa Rao, B. L. S. (1987): Asymptotic Theory of Statistical Inference. Wiley, New York.
[3] Hurt, J. (1986): Asymptotic expansions for moments of functions of stochastic processes and their applications. Statistics & Decisions, Vol. 4, pp. 251-271.
[4] Pugh, E. L. (1963): The best estimate of reliability in the exponential case. Oper. Res. 11/17, pp. 57-61.
[5] Hodges, J. L., jr. and Lehmann, E. L. (1970): Deficiency. Ann. Math. Statist. Vol. 41, pp. 783-801.
[6] Hurt, J. and Wertz, W. (1983): Asymptotic properties of the invariantly optimal estimator of reliability in the 
exponential case. Statistics & Decisions, Vol. 1, pp. 197-204.
[7] Friesl, M. and Hurt, J.  (2007): On Bayesian estimation in an exponential distribution under random censorship. Kybernetika, Vol 43, pp. 45-60.