2
$\begingroup$

FindDistributionParameters is super useful for extracting parameters and using them for other things, like defining random variables.

However I can't find any literature on how to extract the errors in the parameters that FindDistributionParameters determines.

For example:

TestData = RandomVariate[RayleighDistributiom[1.2234],10000]
FindDistributionParameters[TestData, RayleighDistribution[\[Sigma]], ParameterEstimator -> "MaximumLikelihood"]

Now obviously I can just run this $N$ times to get an idea on the error. Or if I have real data where I only have one set, I can bootstrap it and again run $N$ times. But I'd like to know if there is an inbuilt estimation on the parameter errors already?

$\endgroup$
7
  • $\begingroup$ Have you already seen DistributionFitTest[]? $\endgroup$ – J. M.'s ennui Apr 24 '20 at 0:09
  • 1
    $\begingroup$ Depends on how loose you want to be in the definition of "inbuilt". One can use the result from the LogLikelihood function to obtain estimates of the standard errors of the parameters (along with the estimated covariances). $\endgroup$ – JimB Apr 24 '20 at 0:28
  • $\begingroup$ @J.M. I thought this tests only for normally distributed data? Would you be willing to throw together a quick example? $\endgroup$ – Q.P. Apr 24 '20 at 0:28
  • $\begingroup$ @JimB Any method that just uses the initial set of data would be good! I would just prefer to avoid bootstrapping data. $\endgroup$ – Q.P. Apr 24 '20 at 0:29
  • 1
    $\begingroup$ "I thought this tests only for normally distributed data?" - did you look at e.g. the second example under "Basic Examples" in the help file for DistributionFitTest[], where one is testing for a fit with the Pareto distribution? How about the third and succeeding examples under "Scope"? $\endgroup$ – J. M.'s ennui Apr 24 '20 at 0:32
6
$\begingroup$

I'm glad to see you are wanting to obtain estimates of precision (rather than hypothesis testing): "An estimate of a parameter without an associated measure of precision is at best of unknown value." Of course, getting an estimate of standard error or a precise estimate of a parameter doesn't help if the distribution is inappropriate for the data.

The necessary calculations are already performed internally to obtain the maximum likelihood estimate so I don't know why FindDistributionParameters or some other builtin function doesn't offer that. Fortunately the calculations are relatively simple (especially for a one-parameter distribution).

Here's an example that matches your example:

SeedRandom[12345];
TestData = RandomVariate[RayleighDistribution[1.2234], 10000];
mle = FindDistributionParameters[TestData, RayleighDistribution[σ],
  ParameterEstimator -> "MaximumLikelihood"]
(* {σ -> 1.21208} *)

logL = LogLikelihood[RayleighDistribution[σ], TestData];
se = Sqrt[-1/(D[logL, {σ, 2}]) /. mle]
(* 0.0060604 *)

When there are multiple parameters, an Inverse is necessary which results in an estimated covariance matrix. Here's an example for a 2-parameter Weibull:

(* Generate data from a 2-parameter Weibull *)
SeedRandom[12345];
TestData = RandomVariate[WeibullDistribution[1.2345, 5.4321, 10.9876], 10000];

(* Find maximum likelihood estimates *)
mle = FindDistributionParameters[TestData, WeibullDistribution[α, β, 10.9876],
  ParameterEstimator -> "MaximumLikelihood"]
(* {α -> 1.24141, β -> 5.36107} *) 

(* Get the log of the likelihood *)
logL = LogLikelihood[WeibullDistribution[α, β, 10.9876], TestData];

(* Find the hessian evaluated at the maximum likelihood solution *)
hessian = (D[logL, {{α, β}, 2}]) /. mle;

(* Get estimate of covariance matrix *)
(cov = -Inverse[hessian]) // MatrixForm

$$\left( \begin{array}{cc} 0.000093733 & 0.000137995 \\ 0.000137995 & 0.00206814 \\ \end{array} \right)$$

(* Standard errors *)
seα = cov[[1, 1]]^0.5
(* 0.00968158 *)
seβ = cov[[2, 2]]^0.5
(* 0.0454768 *)
$\endgroup$
9
  • $\begingroup$ You're a bloody marvel. Can I ask what the meaning of se = Sqrt[-1/(D[logL, {σ, 2}]) /. mle] is? As in why is the equation of that form, does it have a name so I can read up about it? $\endgroup$ – Q.P. Apr 24 '20 at 1:07
  • $\begingroup$ se represents the estimate of the standard error associated with the parameter estimate. The reciprocal of minus the second derivative of the log of the likelihood evaluated at the maximum likelihood estimate of $\sigma$ obtains the estimate of the variance. To get the standard error one takes the square root of that. However, I suspect you want a general description of the procedure. I'll look for a concise blurb on that. $\endgroup$ – JimB Apr 24 '20 at 1:16
  • $\begingroup$ superb. I shall accept your answer! $\endgroup$ – Q.P. Apr 24 '20 at 1:20
  • $\begingroup$ final question, is the procedure any different for distributions with more than one parameter? I ask as I adapted the code for a three parameter distribution and the derivative hasn't yet compiled: TestData = RandomVariate[WeibullDistribution[1.2345, 5.4321, 10.9876], 10000]; mle = FindDistributionParameters[TestData, WeibullDistribution[\[Alpha], \[Beta], \[Gamma]], ParameterEstimator -> "MaximumLikelihood"] logL = LogLikelihood[WeibullDistribution[\[Alpha], \[Beta], \[Gamma]], TestData]; se = Sqrt[-1/(D[logL, {\[Alpha], 2}, {\[Beta], 2}, {\[Gamma], 2}])] $\endgroup$ – Q.P. Apr 24 '20 at 2:48
  • 1
    $\begingroup$ You've pick an example where no maximum likelihood estimate exists when all 3 of the Weibull parameters are unknown. I'll put in a 2-parameter Weibull example. $\endgroup$ – JimB Apr 24 '20 at 3:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.