# FindDistribution doesn't return a decent result

I have the following data set:

data = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11};

where its histogram is:

Now, I want to find the appropriate underlying distribution, so one can write:

dist = FindDistribution[data]

which returns $$NB(4, 0.4912)$$. However, as the following plot shows, the negative binomial distribution doesn't seem to be a very good choice:

How can I get a more decent distribution?

EDIT

I used SmoothKernelDistribution in my previous post for another seta of data, which gave wrong result. So, I thought it's better to use FindDistribution.

• Nature is under no obligation to operate according to any distribution with convenient mathematical properties. Aug 30 '21 at 17:39

One may consider several variants

data = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7,
7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8,
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9,
10, 10, 10, 10, 11, 11, 11}; FindDistribution[data, 3]

{NegativeBinomialDistribution[4, 0.49117], MixtureDistribution[{0.414276, 0.362911, 0.222813}, {DiscreteUniformDistribution[{0, 3}], DiscreteUniformDistribution[{4, 6}], BinomialDistribution[11, 0.705747]}], MixtureDistribution[{0.414716, 0.38796, 0.197324}, {DiscreteUniformDistribution[{0, 3}], DiscreteUniformDistribution[{4, 6}], PascalDistribution[7, 0.869474]}]}

Show[{Histogram[data, 11, "PDF"],
DiscretePlot[PDF[MixtureDistribution[{0.41427560008093844, 0.362910940968043,
0.22281345895101873}, {DiscreteUniformDistribution[{0, 3}],
DiscreteUniformDistribution[{4, 6}],
BinomialDistribution[11, 0.7057474392420473]}], x], {x, 0,
12}]}]

to this end. Conclusions about the goodness may be drawn by statistical tests, not by plots.

PearsonChiSquareTest[data,NegativeBinomialDistribution[4, 0.49117043121149895]]

6.18699*10^-8

PearsonChiSquareTest[data,MixtureDistribution[{0.41427560008093844, 0.362910940968043, 0.22281345895101873}, {DiscreteUniformDistribution[{0, 3}], DiscreteUniformDistribution[{4, 6}], BinomialDistribution[11, 0.7057474392420473`]}]]

0.0845034

FindDistribution[data, 3, "PearsonChiSquare",  PerformanceGoal -> "Quality"]

{{NegativeBinomialDistribution[4, 0.49117], 6.18699*10^-8}, {MixtureDistribution[{0.414276, 0.362911, 0.222813}, {DiscreteUniformDistribution[{0, 3}], DiscreteUniformDistribution[{4, 6}], BinomialDistribution[11, 0.705747]}], 0.0845034}, {MixtureDistribution[{0.414716, 0.38796, 0.197324}, {DiscreteUniformDistribution[{0, 3}], DiscreteUniformDistribution[{4, 6}], PascalDistribution[7, 0.869474]}], 0.767572}}

shows, MixtureDistribution[{0.414716, 0.38796, 0.197324}, {DiscreteUniformDistribution[{0, 3}], DiscreteUniformDistribution[{4, 6}], PascalDistribution[7, 0.869474]}] is not so bad fit.

• @Leyli: No, it is not so as PearsonChiSquareTest shows. Aug 29 '21 at 11:06
• @Leyli: Try FindDistribution[data, n, prop]. See the documentation to FindDistribution for more info. Aug 29 '21 at 11:26
• @Leyli: Compare these tests on your own. Good luck! Aug 29 '21 at 12:16
• The PearsonChiSquareTest's in your examples use the fitted values of the parameters but that doesn't account for the fine tuning of those parameters to the data.
– JimB
Aug 30 '21 at 16:36

Context matters. To give a better answer one would need to know how the data was collected, a bit of the subject matter, and your objective.

So I am going to assume that you have a random sample from discrete distribution and that you want a parsimonious description of the resulting data.

A complete description of you data consists of the 12 relative frequencies but having 12 parameters is likely not a parsimonious description that would meet your unstated objective.

No common parametric discrete distributions typically produce samples that look like your data (the bimodal shape is the characteristic not inherent in most of those common distributions). The fit found by @user64494 is technically fine but it contains 9 parameters to describe 12 frequencies. That would not be parsimonious.

The large number of zeros suggest that fitting a distribution with an inflated number of zeros might be an option. Again, here is where subject matter context matters. Mathematica can do this with a MixtureDistribution.

Consider a mixture distribution where 100w% of the time one obtains a zero and otherwise one obtains a sample from a negative binomial distribution. The estimates of the parameters proceeds as follows:

d = MixtureDistribution[{w, 1 - w}, {BinomialDistribution[1, 0],
NegativeBinomialDistribution[n, p]}];
sol = FindDistributionParameters[data, d]
(* {w -> 0.104849, n -> 13.6055, p -> 0.748738} *)

nparms = 3; (* Number of parameters estimated *)
aic = -2 LogLikelihood[d /. sol, data] + 2*nparms
(* 1388.29 *)

(* Pearson chisquare test *)
PearsonChiSquareTest[data, d]
(* 0.135872 *)

(* Estimated probability mass function  *)
pmf = Table[{x, PDF[d /. sol, x]}, {x, 0, 11}];

(* Relative frequencies  of data *)
f = Tally[data]
f[[All, 2]] = f[[All, 2]]/Length[data]

(* Estimated pmf and relative frequencies *)
ListPlot[{f, pdf}, PlotLegends -> {"Observed", "Fitted"}]]

(The aic value will be considered later.)

Another alternative is a mixture of two Poisson distributions. This also results in the estimation of just 3 parameters:

d = MixtureDistribution[{w, 1 - w}, {PoissonDistribution[λ1],
PoissonDistribution[λ2]}];
sol = FindDistributionParameters[data, d]
(* {w -> 0.820217, \[Lambda]1 -> 4.93884, \[Lambda]2 -> 0.46591} *)
nparms = 3;
aic = -2 LogLikelihood[d /. sol, data] + 2*nparms
(* 1382.35 *)
PearsonChiSquareTest[data, d]
(* 0.524542 *)
pmf = Table[{x, PDF[d /. sol, x]}, {x, 0, 11}];
f = Tally[data]
f[[All, 2]] = f[[All, 2]]/Length[data]
ListPlot[{f, pmf}, PlotLegends -> {"Observed", "Fitted"}]

Both fits do not suggest a statistically significant departure from the two mixture distributions but the aic is considerably less for the mixture of 2 Poisson distributions.

To summarize my sermons: (1) Subject context matters, (2) Maybe consider just presenting the 12 relative frequencies if there's no justification for attempting to find either a parsimonious description or if there's no theoretical reason that there should be some nice distribution that generates such data, (3) Such questions might be better asked at CrossValidated for statistical issues and then back here for implementation with Mathematica.

• Indeed, all that is not so simple. Here is an example: data = RandomVariate[WeibullDistribution[2, 1.9], 400];FindDistribution[data, 4, "PearsonChiSquare"] results in {{RayleighDistribution[1.28483], 0.676653}, {WeibullDistribution[1.89799, 1.82689], 0.690301}, {WeibullDistribution[1.83899, 1.78166, 0.0367204], 0.662881}, {GammaDistribution[2.84702, 0.569162], 0.189092}}. Aug 31 '21 at 8:10