Lesson #8: Statistical Inference

In this lesson, we will almost exclusivly discuss estimation of a parameter in a parametric model.

Let (Subscript[X, 1]),(Subscript[X, 2]),...,(Subscript[X, n]) be a "random i.i.d. sample" from distribution f((Subscript[x, i])|θ).

We are interested in using the observations (Subscript[X, 1]),(Subscript[X, 2]),...,(Subscript[X, n]) to "guess" what is g(θ), often just θ.

Example:
Let (Subscript[X, i])~Gamma(α,λ).
Here θ = {α,λ}.  We may be interested in estimating g(θ)=(α/λ), the mean. Or,  g(θ)={θ}. Or g(θ)=(α/(λ^2)), the variance, etc..

An estimator, (OverBar[θ]) (I wanted hat)   is a function of (Subscript[X, 1]),(Subscript[X, 2]),...,(Subscript[X, n]).  It is thus a random variable.

Attributes of Estimators

Bias

myVariance[xx_]:=((Table[(((xx[[i]]-Mean[xx]))^2),{i,1,Length[xx]}]//Total)/(Length[xx]))

nnn=10000;

rvs=RandomReal[NormalDistribution[],nnn];

myVariance[rvs]

Table[myVariance[RandomReal[NormalDistribution[],5]],{10000}] //Mean

((5-1)/5)//N

So this is why we use (1/(n-1)) In denominator of sample variance

Consistency

xx=RandomReal[NormalDistribution[],10000];

Table[myVariance[Take[xx,n]],{n,10,1000,10}];

ListPlot[%]

Maximum Likelihood Estimation

Example - Bernoullli(p)

PDF[BernoulliDistribution[p],x]

like[p_,xx_]:=Apply[Times,
Table[(p^(xx[[i]]))(((1-p))^(1-xx[[i]])),{i,1,Length[xx]}]
]

like[p,{1,0,1}]

Plot[%,{p,0,1}]

y is total of xx.  n is the length

logLike[p_,y_]:=y Log[p]+(n-y)Log[1-p]

Clear[xx]

Solve[D[logLike[p,y],p]==0,p]

Tuples[{a,b,c},2]

res=Tuples[{0,1},5]

res[[15]]

posLike=Map[
like[p,#]&,
res]

Manipulate[
Plot[posLike[[i]],{p,0,1}],
{i,1,Length[res],1},SaveDefinitions->True]

Continue to make a demonstration ....

Example - Gamma(α,β)

PDF[GammaDistribution[α,β],x]

(((E^(-(x/β))) (x^(-1+α)) (β^(-α)))/(Gamma[α]))

Mean[GammaDistribution[α,β]]

like[α_,β_,xx_]:=Apply[Times,Table[(((E^(-((xx[[i]])/β))) ((xx[[i]])^(-1+α)) (β^(-α)))/(Gamma[α])),{i,1,Length[xx]}]]

Plot3D[
like[α,β,{2,4,8,2,9}],
{α,0,10},{β,0,5},PlotPoints->100,PlotRange->All,AxesLabel->{α,β}]

plotSample[xx_]:=Plot3D[like[α,β,xx],{α,0,10},{β,0,10},PlotRange->All,AxesLabel->{α,β}]

RandomReal[GammaDistribution[3,3],10]

Table[plotSample[RandomReal[GammaDistribution[3,3],10]],{5}]

Table[plotLogSample[RandomReal[GammaDistribution[3,3],10]],{5}]

plotLogSample[xx_]:=Plot3D[Log[like[α,β,xx]],{α,0,10},{β,0,10},PlotRange->All,AxesLabel->{α,β}]

sample=RandomReal[GammaDistribution[2,2],10];
plotSample[sample]
plotLogSample[sample]

?NIntegrate

?NSolve

?NDSolve

?NMaximize

sample={4,5,4,5,6.3423,7};

NMaximize[like[α,β,{2,4,8,2,9}],{α,β}]

NMaximize[Log[like[α,β,{2,4,8,2,9}]],{α,β}]

NMaximize[Log[like[α,β,sample]],{α,β}]

Alternative : method of moments:

Mean[GammaDistribution[α,β]]

Variance[GammaDistribution[α,β]]

Solve[
{Mean[GammaDistribution[α,β]] == xBar,
Variance[GammaDistribution[α,β]] == sSquared},{α,β}]

sample

estMoments[sample_]:={(((Mean[sample])^2)/(Variance[sample])),((Variance[sample])/(Mean[sample]))}

estMoments[sample]

estMoments[RandomReal[GammaDistribution[2,2],100000]]

Now lets make a "serious" demonstration of this ...

Demonstration During Class

logLike[α_,β_,xx_]:=Log[Apply[Times,Table[(((E^(-((xx[[i]])/β))) ((xx[[i]])^(-1+α)) (β^(-α)))/(Gamma[α])),{i,1,Length[xx]}]]]

Clear[sampleSize];

testAlpha=1;testBetta=1;
repeats=30;
Table[
Map[{α,β} /.NMaximize[#,{α,β}][[2]]&,
Table[logLike[α,β,RandomReal[GammaDistribution[testAlpha,testBetta],sampleSize]],
{30}]
],{sampleSize,100,1000,100}]

data=%100;

Length[data]

Manipulate[ListPlot[data[[i]],AxesOrigin->{0,0},PlotRange->{{0,2},{0,2}},Epilog->{Red,PointSize[Large],Point[{1,1}]}]
,
{i,1,10,1}]

ListPlot[%93,AxesOrigin->{0,0},Epilog->{Red,PointSize[Large],Point[{1,1}]}]

In Class Above

testAlpha=1;testBetta=1;
repeats=30;
sampleSize=2000;
samples=Table[RandomReal[GammaDistribution[testAlpha,testBetta],sampleSize],{repeats}];
logLikes=Map[Log[like[α,β,#]]&,samples];
estsMLE=Map[({α,β}/.NMaximize[#,{α,β}][[2]])&,logLikes];
estsMoment=Map[estMoments,samples];

ListPlot[{estsMLE,estsMoment},AxesOrigin->{0.8,0.8}]