Standard errors for maximum likelihood estimates in FindDistributionParameters

Question

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?

Answer 1

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} *)

Answer 2

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}} *)