3
$\begingroup$

The first line generates 1050 data points to fit. The second line fits this data to the same distribution. This takes about 100ms on my computer. Doing this with FindDistributionParameters does not speed up the fit. Surprisingly, giving initial guesses that match the parent distribution does not speed up the fit. Is there any other way to speed up this code? Can it be compiled?

fakedata = RandomVariate[ExtremeValueDistribution[0.15, 0.11], 1050];
AbsoluteTiming[
 EstimatedDistribution[fakedata, 
   ExtremeValueDistribution[\[Alpha], \[Beta]]];]
$\endgroup$
  • 3
    $\begingroup$ You can try a different ParameterEstimator (and different Methods to these). ParameterEstimator -> "MethodOf*Moments" is incredibly fast, but I assume potentially inaccurate. It appears the max-likelihood estimator is the only slow one (but also probably the best one for real data). $\endgroup$ – b3m2a1 Jul 5 '18 at 23:31
  • $\begingroup$ Thanks, I tried your suggestion and it gives of speedup of about 100x! But I will need to check out the accuracy. $\endgroup$ – Michael B. Heaney Jul 5 '18 at 23:34
  • $\begingroup$ I checked out the accuracy, and it is the same (within the statistical error) for both! $\endgroup$ – Michael B. Heaney Jul 6 '18 at 16:34
  • 1
    $\begingroup$ Just a note stating that my answer does a pretty marginal improvement in speed (like less than 5% faster). One of the reasons is that the maximum likelihood estimators require the ratio of $\sum_{i=1}^n x_i e^{-x_i/\beta}$ and $\sum_{i=1}^n e^{-x_i/\beta}$ which are expensive calculations if needed to be performed for many iterations. See stats.stackexchange.com/questions/71197/…. $\endgroup$ – JimB Aug 14 '18 at 6:33
  • $\begingroup$ Thanks, that is helpful! $\endgroup$ – Michael B. Heaney Aug 14 '18 at 14:50
3
$\begingroup$

Please bear with me as I go through a few steps. First generate some data:

SeedRandom[12345];
fakedata = RandomVariate[ExtremeValueDistribution[0.15, 0.11], 1050];

Method of moments is pretty quick:

AbsoluteTiming[
 mom = EstimatedDistribution[fakedata, 
   ExtremeValueDistribution[α, β],
   ParameterEstimator -> "MethodOfMoments"]]
(* {0.001024, ExtremeValueDistribution[0.150455, 0.105036]} *)

Maximum likelihood takes much longer - as you noticed:

AbsoluteTiming[
 mle = EstimatedDistribution[fakedata, ExtremeValueDistribution[α, β]]]
(*{0.111683, ExtremeValueDistribution[0.150217, 0.106278]} *)

But what if we plug in the method of moments estimators as starting values? That speeds things up a little bit:

AbsoluteTiming[
 EstimatedDistribution[fakedata, ExtremeValueDistribution[α, β],
  {{α, mom[[1]]}, {β, mom[[2]]}}]]
(* {0.0896438, ExtremeValueDistribution[0.150217, 0.106278]} *)

But you really always need estimates of the standard errors which EstimatedDistribution doesn't give you but are relatively easy to get:

logL = LogLikelihood[ExtremeValueDistribution[α, β], fakedata];
cov = -Inverse[(D[logL, {{α, β}, 2}]) /. {α -> mle[[1]], β -> mle[[2]]}];
seα = cov[[1, 1]]^0.5
seβ = cov[[2, 2]]^0.5
$\endgroup$

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.