BeginPackage["MaximumLikelihood`"]; ML::usage= "ML[data,f,startingvalues] returns the maximum likelihood estimates and their \ estimated variance for the probability density f and corresponding data. The \ probability density is a pure function that is mapped over the data."; PlotDensities::usage= "PlotDensities[data,f,g] plots the theoretical (f) and estimated densities \ (g) for a continuous distribution."; ListPlotDensities::usage= "ListPlotDensities[data,f,g] plots the theoretical (f) and estimated \ densities (g) for a discrete distribution."; Plot2DConfidenceRegion::usage= "Plot2DConfidenceRegion[estimate] gives a two-dimensional contour plot of the \ 99%, 95%, and 90% confidence regions centered on the maximum likelihood \ estimates."; ML::nodata="Empty data list in ML routine."; ML::failure="Failure to find maximum of likelihood in ML routine."; Begin["`Private`"]; ML[data_List,f_,startingvalues_List]:= Module[{L,results,parameters,h,vcovresults,SE}, If[Length[data]==0,Message[ML::nodata];Return[];]; L=Apply[Plus,Log[Map[f,data]]]; results=Apply[FindMinimum, FlattenAt[{-L,startingvalues,MaxIterations->500},{2}]]; parameters=Transpose[startingvalues][[1]]; If[Not[Apply[And,Map[NumberQ,parameters/.results[[2]]]]], Message[ML::failure];Return[]]; h=Transpose[(Map[D[Log[Map[f,data]],#]&,parameters])/.results[[2]]]; vcovresults=Inverse[Transpose[h].h]; SE=Table[Sqrt[vcovresults[[i,i]] ],{i,Length[parameters]}]; Print["-----Maximum Likelihood with Gradient Variance Computation-----"]; Print["Number of Observations: "<>ToString[Length[data]] ]; Print["Log(L): "<>ToString[-results[[1]] ] ]; Print[NumberForm[ TableForm[Transpose[{parameters, SE, parameters/SE}]/.results[[2]], TableHeadings->{ ToString[#]& /@ parameters,{"Parameter","Standard Error", "z-statistic"}} ],4]]; Flatten[{LogLikelihood->-results[[1]],results[[2]],vcov->vcovresults}] ]; Plot2DConfidenceRegion[estimate_]:= Module[{parameters,coeff,var,d1,d2,t11,t12,t21,t22,g}, parameters=Drop[Drop[Map[#[[1]]&,estimate],1],-1]; coeff=parameters/.estimate; var=vcov/.estimate; g=PDF[MultinormalDistribution[coeff,var],#]&; d1=Sqrt[var[[1,1]]]; d2=Sqrt[var[[2,2]]]; {t11,t12,t21,t22}= {coeff[[1]]-3*d1,coeff[[1]]+3*d1,coeff[[2]]+ 3*d2,coeff[[2]]- 3*d2}; ContourPlot[g[{t1,t2}],{t1,t11,t12},{t2,t21,t22}, PlotLabel->"99%, 95%, and 90% Confidence Regions", FrameLabel->Map[ToString,parameters], ContourShading->False,PlotPoints->40, Contours->{0.01,0.05,0.1}*g[coeff], Epilog->{PointSize[0.04],RGBColor[1,0,0],Point[coeff]/.estimate}] ]; PlotDensities[data_List,f_,g_]:= Module[{nBins,dx,position, height,width,xx,yy,min,max}, min= Min[data]; max= Max[data]; nBins = Floor[5+(Length[data])^.33*1.5]; dx = (max - min)/nBins; position =Table[i,{i,min+dx/2,max-dx/2,dx}]; height =BinCounts[data, {min, max-.00000001, dx}]/(Length[data]*dx); width=Table[dx,{i,nBins}]; xx= GeneralizedBarChart[ Transpose[ {position, height, width}], BarStyle ->(GrayLevel[1]&), DisplayFunction->Identity]; yy= Plot[{f[t],g[t]},{t,Min[data],Max[data]}, Axes->{True,False}, PlotStyle->{{Thickness[.005]}, {Thickness[.005],Dashing[{.02}]}}, DisplayFunction->Identity]; Show[xx,yy,DisplayFunction->$DisplayFunction]; ]; ListPlotDensities[data_List,f_,g_]:= Module[{xx,yy}, xx:= BarChart[Map[{#[[1]]/Length[data],#[[2]]}&,Frequencies[data]], BarStyle ->(GrayLevel[1]&),DisplayFunction->Identity ]; yy= MultipleListPlot[ Table[f[i],{i,Min[data],Max[data]}], Table[g[i],{i,Min[data],Max[data]}], PlotJoined->{True,True}, Axes->{True,False}, DisplayFunction->Identity, PlotStyle->{ {Thickness[.007]}, {Thickness[.007],Dashing[{.02}]}},DisplayFunction->Identity ]; Show[xx,yy,DisplayFunction->$DisplayFunction]; ]; End[]; EndPackage[];