FindDistribution[] with error bars? (And a ChiSquared Question)

I've got some data defined as follows:

geigerData={466,471,500,442,469,448,481,470,486,472,446,487,454,485,453,447,488,454,523,478,461,473,472,466,458,463,451,455,465,488,488,413,448,453,438,453,452,505,471,467,471,469,508,469,487,472,430,479,473,475,481,468,482,459,460,478,478,472,447,471,469,431,444,485,453,480,469,462,460,485,497,493,524,491,447,483,478,441,455,475,460,429,468,471,468,479,503,482,495,463,479,504,463,459,492,482,459,476,506,464,453,467,491,448,458,444,456,429,477,473,477,451,483,468,463,464,442,469,449,460,455,466,492,493,462,478,461,529,462,461,461,457,460,494,473,503,462,475,487,467,478,455,449,462,445,462,461,456,441,458,450,465,452,464,474,460,478,439,447,428,468,476,445,476,439,533,481,436,476,477,491,462,454,472,484,465,459,497,485,467,457,462,481,475,489,466,462,455,440,446,479,444,480,465,455,486,506,471,494,429}

I am trying to fit Gaussian (Title Normal in Mathematica) and Poisson fits to this data. FindDistribution[] appears to be what I'm looking for:

FindDistribution[geigerData,1,All,TargetFunctions->{NormalDistribution}]

This gives me 0.96 for ChiSquare, which is one of two values I'm looking for to talk about the goodness-of-fit. However, I have two questions:

1. The Documentation states that "PearsonChiSquare" gives "PearsonChiSquareTest p-value." Does that mean that what it gives me is the value for ChiSquare or is it the corresponding p-value? I need both, but I can't tell what's been given to me.
2. My data is supposed to have error bars. If I was making a list plot, this would be as easy as:

x_i +- Sqrt(x_i)

Where the Sqrt(x_i) is the margin of error for the data given by a derivation that's not the focus of this post. I would iterate over the list to create a new one with this equation for each point. But when I do this, while it's supported by ListPlot, FindDistribution gives me an error:

Argument (insert list data here) at position 1 does not have the right format. Data should be a numerical array of depth.

Is there another way I can define a list with error that will work with this equation? NonLinearFit has Error -> that can be used, but I don't see that option here.

All help is appreciated, thank you!

• Have you seen this post? community.wolfram.com/groups/-/m/t/2020806. Also, you want to fit a probability distribution from a random sample from that distribution. That makes NonlinearModelFit inappropriate for your objective.
– JimB
Jul 7, 2020 at 5:47
• Yes. This is why I am trying to use FindDistribution[], as stated in my post. I am trying to get it to work with data weighted with error. This is something I know for certain that NonlinearModelFit does, but I can't figure it out with FinDistribution[]. This is why I made my post. Also, the post you linked to was also me. Nobody had responded to it so I tried my luck elsewhere. Jul 7, 2020 at 6:10
• There's no need to "weight" the data. (Are you getting that idea from a class instructor?) It is true that if you have a single observation AND assume you have a Poisson distribution, then the square root of the count can estimate the standard deviation. But that doesn't mean your count has measurement error. Are you expecting more than Poisson variability where there really is an error in the count?
– JimB
Jul 7, 2020 at 6:17
• Yes. We recorded this data in class, there is going to be error in the number of counts. It increases with N, being Sqrt(N). It's not just an idea from the instructor; it's a requirement. Jul 7, 2020 at 6:27
• Unless, of course, I am just extremely confused. Honestly, I've been at this assignment for so long I wouldn't be surprised if I missed something huge. Jul 7, 2020 at 6:30

To obtain what Mathematica uses as a Chisquare test for testing if data comes from a Poisson distribution, then the following will work for your data?

geigerData = {466, 471, 500, 442, 469, 448, 481, 470, 486, 472, 446,
487, 454, 485, 453, 447, 488, 454, 523, 478, 461, 473, 472, 466,
458, 463, 451, 455, 465, 488, 488, 413, 448, 453, 438, 453, 452,
505, 471, 467, 471, 469, 508, 469, 487, 472, 430, 479, 473, 475,
481, 468, 482, 459, 460, 478, 478, 472, 447, 471, 469, 431, 444,
485, 453, 480, 469, 462, 460, 485, 497, 493, 524, 491, 447, 483,
478, 441, 455, 475, 460, 429, 468, 471, 468, 479, 503, 482, 495,
463, 479, 504, 463, 459, 492, 482, 459, 476, 506, 464, 453, 467,
491, 448, 458, 444, 456, 429, 477, 473, 477, 451, 483, 468, 463,
464, 442, 469, 449, 460, 455, 466, 492, 493, 462, 478, 461, 529,
462, 461, 461, 457, 460, 494, 473, 503, 462, 475, 487, 467, 478,
455, 449, 462, 445, 462, 461, 456, 441, 458, 450, 465, 452, 464,
474, 460, 478, 439, 447, 428, 468, 476, 445, 476, 439, 533, 481,
436, 476, 477, 491, 462, 454, 472, 484, 465, 459, 497, 485, 467,
457, 462, 481, 475, 489, 466, 462, 455, 440, 446, 479, 444, 480,
465, 455, 486, 506, 471, 494, 429};

PearsonChiSquareTest[geigerData, PoissonDistribution[λ], "TestDataTable"]


Here is how to duplicate (essentially duplicate rather than exactly duplicate) what Mathematica does using brute force:

(* Get expected counts with approximately equal proportions in each bin *)
nBins = 17
bins = {-1, Table[InverseCDF[PoissonDistribution[mean], i/17], {i, 16}], ∞} // Flatten
(* {-1, 435, 443, 448, 453, 456, 460, 463, 466, 470, 473, 476, 480, 484, 488, 494, 502, ∞} *)
n = Length[geigerData]
(* 200 *)
expected = n*Table[CDF[PoissonDistribution[mean], bins[[i + 1]]] -
CDF[PoissonDistribution[mean], bins[[i]]], {i, nBins}]
(* {12.8004, 12.4428, 11.0582, 13.623, 9.29103, 13.4504, 10.6759, 10.9689,
14.7215, 10.8412, 10.4493, 13.0259, 11.6982, 10.1644, 12.2187, 11.0108,
11.5594} *)

(* Observed counts *)
observed = Table[Length[Select[geigerData, bins[[i]] < # <= bins[[i + 1]] &]], {i, 17}]
(* {7, 9, 14, 12, 11, 15, 18, 11, 16, 15, 9, 16, 10, 12, 10, 4, 11} *)

(* Chisquare statistic *)
chisq = Total[(observed - expected)^2/expected]
(* 18.1329 *)

(* P-value *)
pvalue = 1 - CDF[ChiSquareDistribution[17 - 2], chisq]
(* 0.255744 *)


To test to see if the data could have been generated by a Poisson distribution (with unknown mean), there is absolutely no need to use the square root of the count as a weight. If you have a reference to such a procedure, I'd really like to see it.

It seems like you need to examine different numbers of bins. Here's the above code modified to look at 4 to 20 bins testing for a normal distribution with your data. (This is not how one determines an optimal number of bins. You'd first need to state what is meant by "optimal".)

results = ConstantArray[{0, 0, 0}, 17];
Do[
mean = Mean[geigerData] // N;
stdev = StandardDeviation[geigerData] // N;
bins = {-∞, Table[InverseCDF[NormalDistribution[mean, stdev], i/nBins],
{i, nBins - 1}], ∞} // Flatten;

n = Length[geigerData];
expected = n*Table[CDF[NormalDistribution[mean, stdev], bins[[i + 1]]] -
CDF[NormalDistribution[mean, stdev], bins[[i]]], {i, nBins}];

(* Observed counts *)
observed = Table[Length[Select[geigerData, bins[[i]] < # <= bins[[i + 1]] &]], {i, nBins}];

(* Chisquare statistic *)
chisq = Total[(observed - expected)^2/expected];

(* P-value *)
pvalue = 1 - CDF[ChiSquareDistribution[nBins - 3], chisq];

results[[nBins - 3]] = {nBins, chisq, pvalue},
{nBins, 4, 20}]

TableForm[results, TableHeadings -> {None, {"# of bins",
"\!$$\*SuperscriptBox[\(χ$$, $$2$$]\) statistic",
"P-value"}}]


• Do you know if it's possible to set the bin count of PearsonChiSquaredTest[]? Jul 8, 2020 at 1:23
• I'm not aware of anyway to change the bin count for PearsonChiSquareTest. If it's any consolation, the bin count increases with sample size (as it really should) and can be approximated with the formula Round[2.23117 n^0.385798]. But the above code gets the same results and you can pick any number of bins that you want.
– JimB
Jul 8, 2020 at 2:10
• If I were to use the code above for other types of distributions, such as the Normal and Binomial, would I only need to adjust the distribution used in the beginning and then the degrees of freedom at the end? Jul 8, 2020 at 2:49
• My answer would be "Yes" with a slight hesitation. As you've hinted at the degrees of freedom would be the number of bins minus 1 for matching the total and 1 for every parameter estimated. The "slight hesitation" is that the example data clearly isn't normal as the data is all integers. My professional bias/opinion is testing for any distribution is not very useful. One is usually looking either for a reasonable approximation (estimation rather than hypothesis testing) or looking for gross departures from a distribution.
– JimB
Jul 8, 2020 at 3:28
• Would this work if both expected and actual were both in probability form? Jul 8, 2020 at 3:43