Problem 3 Quiz 7 part(a)
by Nasser Abbasi

Part A

Open the gamma file and load the data

Remove["Global`*"]

Decide on number of bins and make a histogram (see appendix for function I wrote to make probability histogram)

nBins=100;
gz=nmaMakeDensityHistogram[data,nBins];
histPlot=GeneralizedBarChart[gz,BarStyle->White,ImageSize->450,PlotRange->{{Min[data],Max[data]},All},PlotLabel->"Histogram of gamma arrival time data"]

Let me plot a Gamma distribution with the mean arrival time of the data

β=Mean[data]

α=1;

Plot[PDF[GammaDistribution[α,β],x],{x,0,Max[data]},PlotRange->All,PlotLabel->"Gamma[α,β]"]

From the above we see that a Gamma distribution is possible.

Part B

Start by finding the first and second moments of the Gamma distribution

Clear[x,α,β]
m1=ExpectedValue[x,GammaDistribution[α,β],x]

m2=ExpectedValue[(x^2),GammaDistribution[α,β],x]

Now find the first and second moments of the data itself

m1ForData=Mean[data]

m2ForData=Variance[data]+(((Mean[data]))^2)

Now solve for α, β

sol=First@Solve[{m1==m1ForData,m2==m2ForData},{α,β}];
αByMoments=α/.sol

βByMoments=β/.sol

Hence we find λ

Now find estimates of the parameters by MLE method and do the same as above

Clear[α,β]
n=Length[data]

xBar=Mean[data]

eq1=n Log[α] - n Log[xBar]-n PolyGamma[0,α]+(∑_(i=1)^n)Log[data[[i]]]

Solve the above numerically using FindRoot

FindRoot[eq1,{α,1}];
αMLE=α/.%

Now that we found MLE for alpha, use it to find MLE for λ

λMLE=((αMLE)/xBar)

βMLE=(1/(λMLE))

Conclusion for part B

Make a table to compare the α, λ found by the above 2 methods

TableForm[{{αByMoments,"\t",βByMoments},{αMLE,"\t",βMLE}},TableHeadings->{{"Method of moments","MLE"},{"α","","β"}}]

Part C

Put the 2 fitting on top of each others to compare

pdfByMomentMethod=Plot[PDF[GammaDistribution[αByMoments,βByMoments],x],{x,0,Max[data]},PlotRange->All,PlotLabel->"Gamma[α,β]",PlotStyle->{Red,Thick}];
p1=Show[{histPlot,pdfByMomentMethod},PlotLabel->"Fitting Gamma[α,β] by method of moments on data",ImageSize->300];

pdfByMLE=Plot[PDF[GammaDistribution[αMLE,βMLE],x],{x,0,Max[data]},PlotRange->All,PlotLabel->"Gamma[α,β]",PlotStyle->{Red,Thick}];
p2=Show[{histPlot,pdfByMLE},PlotLabel->"Fitting Gamma[α,β] by method of MLE",ImageSize->300];
Grid[ {{p2,p1}}]

Conclusion for part (C)

This plot shows more closely the fitting on top of each others.They are very close so hard to see the difference other than near the high frequency part.

pdfByMomentMethod=Plot[PDF[GammaDistribution[αByMoments,βByMoments],x],{x,0,Max[data]},PlotRange->{{0,300},All},PlotLabel->"Gamma[α,β]",PlotStyle->{Dashed,Thick}];
pdfByMLE=Plot[PDF[GammaDistribution[αMLE,1/λMLE],x],{x,0,Max[data]},PlotRange->{{0,300},All},PlotLabel->"Gamma[α,λ]",PlotStyle->{Red,Thick}];
p2=Show[{histPlot,pdfByMomentMethod,pdfByMLE},PlotLabel->"Fitting Gamma[α,β] by method of MLE and Moments",ImageSize->300,PlotRange->{{0,200},All}]

Part (d)

I use bootstrap to compare the standard error of parameters estimated using both the MLE and Method of Moments.
In Bootstrap, we find the estimate of the parameters from the data itself. Then we generate random data from the same distribution of the same size as the data itself using these original estimates, then we estimate the parameters again from the generated random data. We repeat this process keeping track each time of the parameters estimated from each sample. At the end we find the variance and the mean of these parameters and also make histogram to compare the standard errors from the MLE and the moments method. The method which has smaller standard error will be better.

Write 2 functions which find the parameters estimates using the MLE and the moments method

Remove["Global`*"]

m1=ExpectedValue[x,GammaDistribution[α,β],x];
m2=ExpectedValue[(x^2),GammaDistribution[α,β],x];

getMethodOfMomentsParameters[data_]:=Module[{sol,αByMoments,βByMoments},
m1ForData=Mean[data];
m2ForData=Variance[data]+(((m1ForData))^2);
sol=First@Solve[{m1==m1ForData,m2==m2ForData},{α,β}];
αByMoments=α/.sol;
βByMoments=β/.sol;

{αByMoments,βByMoments}
]

getMLEParameters[data_]:=Module[{sol,xBar,αMLE,λMLE,eq,α,n},
xBar=Mean[data];
n=Length[data];
eq=n Log[α] - n Log[xBar]-n PolyGamma[0,α]+(∑_(i=1)^n)Log[data[[i]]];
sol=FindRoot[eq,{α,1}];
αMLE=α/.sol;
λMLE=((αMLE)/xBar);

{αMLE,1/λMLE}
]

getStandardError[numberOfIteration_,data_,sizeOfData_,f_]:=Module[{},
alpha=Table[0,{i,numberOfIteration}];
beta=Table[0,{i,numberOfIteration}];

{dataAlpha,dataBeta}=f[data];

For[i=1,i≤numberOfIteration,i++,
{
sample=RandomReal[GammaDistribution[dataAlpha,dataBeta],sizeOfData];
{alpha[[i]],beta[[i]]}=f[sample];
}
];
{Sqrt[Variance[alpha]],Mean[alpha],Sqrt[Variance[beta]],Mean[beta]}
]

Now call the above functions to find the standard errors

SeedRandom[01010101];

file="E:/nabbasi/data/nabbasi_web_Page/my_courses/FULLERTON_COURSES/Fall_2007/math_502_probability_and_statistics/quiz/quiz7/gamma-arrivals.txt";

data=Flatten[Import[file,"Table"]];
sizeOfData=Length[data]; (*data size*)
numberOfIterations=5000;

resultMoments=getStandardError[numberOfIterations,data,sizeOfData,getMethodOfMomentsParameters]

Now do the same for MLE method

resultMLE=getStandardError[numberOfIterations,data,sizeOfData,getMLEParameters];

Summary table for part D

TableForm[{  {resultMoments[[1]]," ",resultMoments[[2]]," ",resultMoments[[3]]," ",resultMoments[[4]]},
{resultMLE[[1]]," ",resultMLE[[2]]," ",resultMLE[[3]]," ",resultMLE[[4]]}},TableHeadings->{{"Method of moments","MLE"},{"std(α)"," ","Mean(α)"," ","std(β)"," ","Mean(β)"}}]

Conclusion for part D

We see from the table above that standard error for the MLE method for the parameter α is a little smaller than that for the moments method. Also for the β parameter, the MLE method gives a smaller standard error than the moments method.

We conclude that the MLE is better than the moments method since its standard error is smaller. The difference is not too large here in standard error, I was expecting a larger error than this. May be if we have larger data size we can see this.

Part E

Now we need to find confidence intervals for both α and β based on the MLE and the moments methods estimates. Hence we will need to generate 4 C.I's since we have 2 parameters and 2 methods. I will use the 95% C.I. in all cases.
This C.I. is then a random variable which is said to contain the parameter being estimated with probability .95
There are 3 methods to determine the CI for the method of moments and the MLE. Exact methods, approximation based on large sample size, and bootstrap methods. We are asked to use the bootstrap method here.

Algorithm

We start first by writing down the algorithm for finding CI (95 %) for MLE or Moments methods. This algorithm can ofcourse be used for other CI other than 95 % by changing the quantile value, but I used 95 % here for illustration only

1) Load the data from file
2) use the method of MLE to estimate the parameters. Let the estimated parameters be (OverHat[α]) and (OverHat[β])
3) Now, using (OverHat[α])and(OverHat[β]) generate 1000 random samples each of size equal to the original data. Use the distribution Gamma.
4) From each sample generated in step (3) determine using MLE an estimate of the sample parameters α,β call these (SuperStar[(Subscript[α, j])]) and (SuperStar[(Subscript[β, j])]) where j is the sample number
5) Sort the parameters sequence (SuperStar[α]) and the parameters sequence (SuperStar[β]) from small to large
6) Find the parameter (SuperStar[α]) at position 25 and at position 975. Call these (Subscript[(SuperStar[α]), 25]) and (Subscript[(SuperStar[α]), 975]) and do the same for (SuperStar[β]) and Call these (Subscript[(SuperStar[β]), 25]) and (Subscript[(SuperStar[β]), 975])
7)Let (UnderBar[δ])=(Subscript[(SuperStar[α]), 25])-(OverHat[α]) and let (Overscript[δ, -])=(Subscript[(SuperStar[α]), 975])-(OverHat[α]) , Hence the 95% CI for (OverHat[α])=((OverHat[α])-(UnderBar[δ]),(OverHat[α])-(Overscript[δ, -]))
8)Let (UnderBar[δ])=(Subscript[(SuperStar[β]), 25])-(OverHat[α]) and let (Overscript[δ, -])=(Subscript[((SuperStar[β])), 975])-(OverHat[α]) , Hence the 95% CI for (OverHat[β])=((OverHat[β])-(UnderBar[δ]),(OverHat[β])-(Overscript[δ, -]))

Code Implementation

SeedRandom[01010101];

file="E:/nabbasi/data/nabbasi_web_Page/my_courses/FULLERTON_COURSES/Fall_2007/math_502_probability_and_statistics/quiz/quiz7/gamma-arrivals.txt";

data=Flatten[Import[file,"Table"]];
sizeOfData=Length[data]; (*data size*)
numberOfIterations=1000;

getCI[data_,f_]:=Module[{},
alphaStar=Table[0,{i,numberOfIterations}];
betaStar=Table[0,{i,numberOfIterations}];

{alphaHat,betaHat}=f[data];

For[i=1,i≤numberOfIterations,i++,
{
sample=RandomReal[GammaDistribution[alphaHat,betaHat],sizeOfData];
{alphaStar[[i]],betaStar[[i]]}=f[sample];
}
];
alphaStar=Sort[alphaStar];
betaStar=Sort[betaStar];
lowerQAlpha=alphaStar[[25]];
upperQAlpha=alphaStar[[975]];
lowerQBeta=betaStar[[25]];
upperQBeta=betaStar[[975]];

{alphaHat-(upperQAlpha-alphaHat),alphaHat-(lowerQAlpha-alphaHat),betaHat-(upperQBeta-betaHat),betaHat-(lowerQBeta-betaHat)}
]

{alphaLow,alphaHigh,betaLow,betaHigh}=getCI[data,getMethodOfMomentsParameters];

Print["95% C.I. for α using Method of moments is ("<>ToString[alphaLow]<> ","<>ToString[alphaHigh]<>")"];
Print["95% C.I. for β using Method of moments is ("<>ToString[betaLow]<> ","<>ToString[betaHigh]<>")"];

Now do the same for the MLE

{alphaLow,alphaHigh,betaLow,betaHigh}=getCI[data,getMLEParameters];

Print["95% C.I. for α using MLE is ("<>ToString[alphaLow]<> ","<>ToString[alphaHigh]<>")"];
Print["95% C.I. for β using MLE is ("<>ToString[betaLow]<> ","<>ToString[betaHigh]<>")"];

Conclusion for part E

The confidence interval using MLE for the same percentage is smaller than that for the method of moments. This is for both the α and β.
This is good. The smaller the CI is for the same confidence the better it is.  This shows that the MLE method of estimation is better than the method of moments (for large sample size as the case here).

Part f

To answer this part, we use MLE to try to fit Poisson distribution on the data. We use MLE to estimate the λ parameter for the poisson distribution and plot the pdf over the data histogram to see if the fit is good
We know that the MLE estimate of λ for a Poisson is (OverBar[X])

Remove["Global`*"]

lambdaHat=Mean[data];
dist=PDF[PoissonDistribution[lambdaHat],t];
poissonData=Table[{data[[i]],dist/.t->data[[i]]},{i,1,Length[data]}];
p1=ListPlot[poissonData,Joined->False,PlotRange->All,PlotLabel->"Fitting Poisson on the data"]

Lets superimpose the histogram

nBins=100;
gz=nmaMakeDensityHistogram[data,nBins];
histPlot=GeneralizedBarChart[gz,BarStyle->White,ImageSize->450,PlotRange->{{Min[data],Max[data]},All},PlotLabel->"Histogram of gamma arrival time data"];

Show[{histPlot,p1},PlotRange->All]

Conclusion part f

From the above we see that a Poisson distribution is not a good probability law to describe the data. Poisson fits a random variable which represents the number of times an event occur in some fixed time. The parameter λ represents the average rate at which the event occur per unit time.
The interarrivals times (from the histogram of the data itself) shows that more events occur (particles arrive) which has the time between them small than when the time is large between them. But in Poisson, the rate should be fixed and not changing. This explains why Poisson does not give a good fit.
(I really would like to know the correct answer for this one if mine is not correct. Thank you)

Appendix

Contain functions needed

file="E:/nabbasi/data/nabbasi_web_Page/my_courses/FULLERTON_COURSES/Fall_2007/math_502_probability_and_statistics/quiz/quiz7/gamma-arrivals.txt";
data=Flatten[Import[file,"Table"]];

<<"BarCharts`";
nmaMakeDensityHistogram[originalData_,nBins_]:=Module[{freq,binSize,from,to,scaleFactor,j,a,currentArea},to=Max[originalData];from=Min[originalData];
binSize=(to-from)/nBins;freq=BinCounts[originalData,binSize];currentArea=Sum[binSize*freq[[i]],{i,nBins}];
freq=freq/currentArea;a=from;
Table[{a+(j-1)*binSize,freq[[j]],binSize},{j,1,nBins}]];