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?
$\begingroup$
$\endgroup$
6
4 Answers
$\begingroup$
$\endgroup$
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} *)
-
1$\begingroup$ 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! $\endgroup$– matildeJun 6, 2012 at 20:34
$\begingroup$
$\endgroup$
10
A possibility is to use a parametric bootstrap. It works as follows:
- Based on the original dataset an underlying distribution is estimated. We want to know the standard deviation (SD) of the resulting parameters.
- From this parametric distribution we randomly draw a large number of datasets similar to the original dataset.
- 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.
- 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}} *)
answered Jun 6, 2012 at 20:21
-
$\begingroup$ excellent, I will try this too to see if it makes a big difference in my sample $\endgroup$– matildeJun 6, 2012 at 20:37
-
$\begingroup$ @matilde Added confidence intervals too. $\endgroup$ Jun 6, 2012 at 20:41
-
1$\begingroup$ @SjoerdC.deVries Isn't this a good use-case for
ParallelTable? $\endgroup$– Mr AlphaJun 6, 2012 at 21:33 -
$\begingroup$ @MrAlpha It is. On my 4-core PC it cuts the time in half. $\endgroup$ Jun 6, 2012 at 21:37
-
1$\begingroup$ related library.wolfram.com/infocenter/MathSource/4272 $\endgroup$ Jun 6, 2012 at 23:04
$\begingroup$
$\endgroup$
2
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.
-
$\begingroup$ @Szabolcs: Thank you very much for your editing! $\endgroup$ Jun 20, 2014 at 17:35
-
1$\begingroup$ 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. $\endgroup$ Jun 13, 2015 at 11:43
$\begingroup$
$\endgroup$
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} *)
NonlinearModelFitdoes provide those, but it is less suited for distributions. $\endgroup$NonlinearModelFitworksunder the assumption that the original Yi are independent normally distributed with mean Yi and common standard deviation(sic) $\endgroup$