# 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 Jun 6, 2012 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 Jun 6, 2012 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. Jun 6, 2012 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) Jun 6, 2012 at 18:57
• If you dont have well recognised distribution then how we generate random no???What is Bootstraping? and what is Darren SD
– user29784
May 27, 2015 at 9:56

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! Jun 6, 2012 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 *)
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 Jun 6, 2012 at 20:37
• @matilde Added confidence intervals too. Jun 6, 2012 at 20:41
• @SjoerdC.deVries Isn't this a good use-case for ParallelTable? Jun 6, 2012 at 21:33
• @MrAlpha It is. On my 4-core PC it cuts the time in half. Jun 6, 2012 at 21:37
• Jun 6, 2012 at 23:04

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 *)
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! Jun 20, 2014 at 17:35
• This is a nice alternative to my parametric bootstrap. In my work I use both types. I believe each has its own use case. If you're sure about the type of underlying distribution I'd pick the parametric bootstrap. If the fitted function is only a rough guess AND you have sufficient data samples I'd use the one in your answer. Jun 13, 2015 at 11:43

If you have a case where the distributions of the maximum likelihood estimators are not very Gaussian-like, then the answer from @SjoerdC.deVries is the way to go.

Otherwise the following will get you estimates of the variances (and covariances).

(* Define a probability distribution with a pdf *)
d = ProbabilityDistribution[Exp[-(x - μ)^2/(2 σ^2)]/Sqrt[2 π σ^2], {x, -∞, ∞},
Assumptions -> σ > 0];

(* Generate some samples from that distribution with specified parameters *)
SeedRandom[12345];
z = RandomVariate[d /. {μ -> 5, σ -> 3}, 100];

(* Get maximum likelihood estimates of parameters *)
mle = FindDistributionParameters[z, d]
(* {μ -> 4.88313, σ -> 2.51668} *)

(* Construct the log likelihood function *)
logL = LogLikelihood[d, z];

(* Obtain estimated covariance matrix using observed Fisher information *)
parameters = {μ, σ};
observedFisherInformation = -D[logL, {parameters, 2}] /. mle;
(cov = Inverse[observedFisherInformation]) // MatrixForm


$$\left( \begin{array}{cc} 0.0633366 & \text{5.051329868355902\grave{ }*{}^{\wedge}-13} \\ \text{5.051222980001161\grave{ }*{}^{\wedge}-13} & 0.0316683 \\ \end{array} \right)$$

(* Standard errors for μ and σ *)
Sqrt[Diagonal[cov]]
(* {0.251668, 0.177956} *)