0
$\begingroup$

I want to extract the distribution parameters for a set of data. I know what the distribution is. However, I find that for different sets or even subsets of this data the values of the extracted distribution using FindDistributionParameters[Data, Distribution] fluctuate quite a lot.

To test this I made a small simulation:

TestData = RandomVariate[WeibullDistribution[84.2, 366.9, -494.9], 10000];

TestDistParams = 
  {αW, βW, γW} /. 
    FindDistributionParameters[
      TestData, 
      WeibullDistribution[αW, βW, γW], 
      {{αW, 84.2}, {βW, 366.9}, {γW, -494.9}}]

Show[
  {Histogram[TestData, {"Raw", Round[Sqrt[Length[TestData]]]}],
   Plot[
     PDF[WeibullDistribution[TestDistParams[[1]], TestDistParams[[2]], TestDistParams[[3]]]][x], 
     {x, -160, -115}]}]

What I find is that the output of FindDistributionParameters[..] fluctuates quite a lot, sometimes the extracted parameters are twice as large as those input into generating the test distribution -- sometimes more so. This is even case when simulating with a large sample (10000 points) and when initialising FindDistributionParameters with guesses exactly as that used to generate the distribution.

What I have found to be far more robust is actually just finding the values of the centre bins with the corresponding PDF value or count, and using NonlinearModelFit to this.

What is the best approach of accurately extracting distribution parameters from a data set, reliably and accurately? Can one use constraints in a similar way to fitting?

On the advice of an experience user, I am adding a specific example. If I use a fixed set of random numbers with SeedRandom[123456]

Using the same code as above I get for SeedRandom[123456] --> {148.383, 653.1, -781.061}. This one is so bad that the plotted PDF is just completely off.

or if I choose another seed, say SeedRandom[851] --> 5.54576*10^6, 2.41777*10^7, -2.41778*10^7

Both examples are for a sample set of 10000 points.

To rephrase my question a little more carefully

My question specifically relates to dealing with FindDistributionParameters when the results are clearly too far off to be considered as part of normal statistical fluctuation (see above examples with fixed seeds) but when the data itself definitely reflects the distribution one is trying to match to it. I.e. when it is specifically drawn from that distribution. Are there constraints one can use for example?

Improve this question
$\endgroup$
12
  • 1
    $\begingroup$ While Mathematica can be used to address this question, it is more appropriate to ask it at CrossValidated (and not at the Mathematics Stack Exchange). Constructing a histogram followed by a regression loses information and is basically a silly way to go. (But, unfortunately you're not alone.) $\endgroup$
    – JimB
    Commented Oct 28, 2019 at 14:06
  • 5
    $\begingroup$ I'm voting to close this question as off-topic because it belongs on CrossValidated. $\endgroup$
    – JimB
    Commented Oct 28, 2019 at 14:07
  • $\begingroup$ But why would they tell me anything different that isn't already implemented in FinDistributionParameters? Assuming Wolfram knows something of statistics, and likelihood evaluations. I am asking specifically why my current approach/use of Mathematica functions doesn't give sensible results. $\endgroup$
    – user27119
    Commented Oct 28, 2019 at 14:21
  • 1
    $\begingroup$ I have no doubt that Wolfram has statisticians far superior to me. However, to leave out associated estimates of precision from FindDistributionParameters shows that decision was made by someone only vaguely familiar with the subject of Statistics. $\endgroup$
    – JimB
    Commented Oct 28, 2019 at 14:45
  • 2
    $\begingroup$ I've retracted my close vote because I now believe there is a numerical instability issue because of the large values being considered for the $\alpha$ parameter of the Weibull distribution. (I'll work on an example that hopefully supports that numerical instability hypothesis.) And there might still be some too high expectations of how well a large sample will do. $\endgroup$
    – JimB
    Commented Oct 28, 2019 at 18:48

1 Answer 1

Reset to default
1
$\begingroup$

This is an extended comment.

I know you believe that the software is misbehaving. However, your example does not show a poor fit to the density. See below.

SeedRandom[123456];
x = RandomVariate[WeibullDistribution[148.383, 653.1, -781.061], 10000];
sol = FindDistributionParameters[x, WeibullDistribution[α, β, μ]]
(* {α -> 668.203, β -> 2970.61, μ -> -3098.53} *)

Show[Histogram[x, Automatic, "PDF"],
  Plot[{PDF[WeibullDistribution[148.383, 653.1, -781.061], z],
    PDF[WeibullDistribution[α, β, μ] /. sol, z]},
    {z, Min[x], Max[x]}, PlotLegends -> {"True", "Estimated"}]]

Histogram with true and estimated density

Notice that one can't tell the true from the estimated density.

Now, you say, what about the estimates of the parameters being "far" from the true values? (I'll deal with that after breakfast.)

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for the extended comment. I have added a couple of other examples using a specific seed where the values are just way off. Those are the ones I want to deal with i.e. can you use constraints in a sensible way if you know enough about your problem? For the case you show I know that pretty much the only way to deal with such errors is bootstrapping and I'm okay with this. It is issue of when FindDistributionParameters routines fails in such a spectacular way that I want to deal with. $\endgroup$
    – user27119
    Commented Oct 28, 2019 at 14:46
  • $\begingroup$ My belief is that FindDistributionParameters is not failing. It is your expectation of the variability associated with the estimators is what needs fixing. That "re-education" is why I have suggested (many times on this forum) to ask at CrossValidated. $\endgroup$
    – JimB
    Commented Oct 28, 2019 at 14:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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