# Standard errors for maximum likelihood estimates in FindDistributionParameters

I am trying to fit a non-standard PDF to data and FindDistributionParameters works great, and gives me the parameters of the distribution back using maximum likelihood, but I cannot find a way to get Mathematica to spit out the standard errors (or confidence intervals) for those estimated parameters. Anyone knows how to do it, or should I just go to a different program, i.e. MATLAB?

-
See, DistributionFitTest[],StatisticalModelAnalysis, and HypothesisTests in the documentation as a start –  Jagra Jun 6 '12 at 17:53
I did look at those as well, but they give me various types of tests on the fit of the distribution, i.e. whether it is normal, etc. and not standard errors, but thanks for your reply –  matilde Jun 6 '12 at 18:04
I think you're right. Apart from p-values for the fit I don't see a way to get more on the parameters of the fit. NonlinearModelFit does provide those, but it is less suited for distributions. –  Sjoerd C. de Vries Jun 6 '12 at 18:29
@SjoerdC.deVries Also, one should always keep in mind that NonlinearModelFit works under the assumption that the original Yi are independent normally distributed with mean Yi and common standard deviation(sic) –  belisarius Jun 6 '12 at 18:57

Under some conditions which are outlined on wikipedia's maximum likelihood page, Sqrt[n](theta-theta0) is asymptotically multivariate normal with mean vector 0 and covariance based on Fisher information. In the formula, n is sample size, theta is the maximum likelihood estimate for the parameter vector, and theta0 is the true (but unknown to us) value of the parameter.

Here is some code that will compute these asymptotic standard errors (provided the log-likelihood is symbolically differentiable).

dist = GammaDistribution[a, b];
data = RandomVariate[dist /. {a -> 3, b -> 10}, 1000];

est = FindDistributionParameters[data, dist]

(* output: {a -> 3.14198, b -> 10.0998} *)

standarderrors[data_, dist_, paramlist_, mleRule_] :=
Block[{len, infmat, cov},
len = Length[data];
(* compute negative of expected Fisher information *)
infmat = -D[LogLikelihood[dist, data], {paramlist, 2}]/len /. mleRule;
(* invert to get asymptotic covariance for Sqrt[n](theta-theta0) *)
cov = Inverse[infmat];
(* standard errors are the Sqrt of diagonal elements divided by sample size *)
Sqrt[Diagonal[cov]/len]
]

standarderrors[data, dist, {a, b}, est]

(* output: {0.133727, 0.466096} *)

-
Wow, thanks a lot for sending me the code. I was going to write it up myself, so you saved me a lot of time! I am a bit surprised that Mathematica would calculate ML estimates, but not provide standard errors, like nonlinearmodelfit does. In any case thanks again! –  matilde Jun 6 '12 at 20:34

A possibility is to use a parametric bootstrap. It works as follows:

1. Based on the original dataset an underlying distribution is estimated. We want to know the standard deviation (SD) of the resulting parameters.
2. From this parametric distribution we randomly draw a large number of datasets similar to the original dataset.
3. Those datasets are then fitted again to the same parametric model, yielding new parameter estimates. We get a new parameter set for each resampled dataset.
4. The SDs of the set of bootstrapped parameters are then used as the SDs of original parameter estimations.

Using Darren's example:

(* Generate a dataset following a GammaDistribution *)
dist = GammaDistribution[a, b];
data = RandomVariate[dist /. {a -> 3, b -> 10}, 1000];

(* Fit distribution parameters *)
est = FindDistributionParameters[data, dist]

(* ==> {a -> 2.912762177, b -> 10.19252432} *)

(* Generate new datasets using the fitted parameters and fit again *)
bootstrap =
Table[
data = RandomVariate[dist /. est, 1000];
{a, b} /. FindDistributionParameters[data, dist],
{10000}
];

(* Darren's SD estimation *)
standarderrors[data, dist, {a, b}, est]

(*   ==> {0.123532062, 0.4717253733}  *)

(* Bootstrapped SD estimation *)
StandardDeviation /@ Transpose[bootstrap]

(*  ==> {0.1242242286, 0.4726653847}   *)


Pretty close, eh?

The advantage of this method is that the conditions for this to be valid are weaker than Darren's. It is conceptually easier, and requires just two lines of code. The disadvantage is that it is slow (about 20 s in this example, which can be reduced to 10 seconds as Mr.Alpha rightly remarks using ParallelTable instead of Table).

Generating 95% confidence intervals is pretty easy too:

{Quantile[#, .025], Quantile[#, .975]} & /@  Transpose[bootstrap]

(* ==> {{2.688031335, 3.176166312}, {9.277990406, 11.11476555}} *)

-
excellent, I will try this too to see if it makes a big difference in my sample –  matilde Jun 6 '12 at 20:37
@matilde Added confidence intervals too. –  Sjoerd C. de Vries Jun 6 '12 at 20:41
thanks a lot for the help –  matilde Jun 6 '12 at 20:55
@SjoerdC.deVries Isn't this a good use-case for ParallelTable? –  Mr Alpha Jun 6 '12 at 21:33
@MrAlpha It is. On my 4-core PC it cuts the time in half. –  Sjoerd C. de Vries Jun 6 '12 at 21:37

The following code is a minor contribution to the idea given by Mr. Sjoerd C. de Vries. After reading HOW TO|Perform a Bootstrap Analysis in the Mathematica Help Browser, I just found a code to deal with the same thing in a slightly different way:

(* Generate a dataset following a GammaDistribution *)
dist = GammaDistribution[a, b];
data = RandomVariate[dist /. {a -> 3, b -> 10}, 1000];

(* Fit distribution parameters *)
Parametersfit = FindDistributionParameters[data, dist]
(* Output: {a -> 3.13025, b -> 9.35438} *)

(* Bootstrapping *)
Bootstrapping := {a, b} /.
FindDistributionParameters[RandomChoice[data, Length@data], dist];
BootEstimates = ParallelTable[Bootstrapping, {10000}];

(* Bootstrapped SD estimation *)
BootStd = StandardDeviation /@ Transpose[BootEstimates]
(* Output: {0.133276, 0.429775} *)


It took about 11 seconds to complete "Bootstrapping" part.

-
@Szabolcs: Thank you very much for your editing! –  user14070 Jun 20 at 17:35