# Actual bins used for Pearson Chi Square Test

PearsonChiSquareTest[data, dist] computes bins for data as a step towards calculating the test statistic. For continuous distributions, the documentation describes the general approach. If I want to get the actual bins used, particularly for discrete distributions, can it be discovered?

Some code to get started:

SeedRandom;
dist = PoissonDistribution[1.5];
n = 100;
data = RandomVariate[dist, n];
testData = PearsonChiSquareTest[data, dist, "TestData"]
(* {3.31031, 0.507301}  *)
nbins = 2 n^(2/5) // Ceiling;
del = DeleteDuplicates[Quantile[dist, Range[0, 1, 1/nbins]]];
observedFrequency = BinCounts[data, {del}];
expectedFrequency = n*(CDF[dist, Most[del]] - CDF[dist, Most[del - 1]]);
chisq = Total[(observedFrequency - expectedFrequency)^2/expectedFrequency]
chisq == testData[]
(* False *)

• There is no mechanism I'm aware of to pull out the bins. – ciao Jun 2 '15 at 5:18

## StatisticsGoodnessOfFitTestingDumppearsonChiSquareTest

Key calculations for PearsonChiSquareTest are performed inside the function StatisticsGoodnessOfFitTestingDumppearsonChiSquareTest. SeedRandom;
mp = MachinePrecision;
dist = PoissonDistribution[1.5];
n = 100;
data = RandomVariate[dist, n];
testData = PearsonChiSquareTest[data, dist, "TestData"]


{3.31031, 0.507301}

With the input (data, null distribution, precision, "SampleToNull") it returns the Pearson test statistic and the observed counts:

StatisticsGoodnessOfFitTestingDumppearsonChiSquareTest[data, dist,  mp, "SampleToNull"]


{3.31031, {25, 36, 27, 8, 4}}

The following lines are the relevant parts of the code for discrete univariate distributions:

StatisticsLibraryDiscreteUnivariateDistributionQ[dist]


True

Bins

nbins = StatisticsGoodnessOfFitTestingDumpiOptimalBinCount[n]


13

bins = Quiet[N[Sort[Quantile[dist, N[Range[0, 1, 1/nbins], mp]], Less], mp]];
bins = bins + N[1/2, mp];
bins[] = bins[] - 1;
bins


{-0.5, 0.5, 0.5, 1.5, 1.5, 1.5, 1.5, 1.5, 2.5, 2.5, 2.5, 3.5, 3.5, 18.5}

Frequencies:

cdfs = StatisticsGoodnessOfFitTestingDumpiCDF[dist, bins];
cdfs = Rest[cdfs] - Most[cdfs];
exp = N[n cdfs, MachinePrecision]


{22.313, 0., 33.4695, 0., 0., 0., 0., 25.1021, 0., 0., 12.5511, 0., 6.56425}

obs = (Length[Cases[data, x_ /; x < #1]] &) /@ bins;
obs = Rest[obs] - Most[obs]


{25, 0, 36, 0, 0, 0, 0, 27, 0, 0, 8, 0, 4}

locs = Position[exp, 0 | 0.];
exp = Delete[exp, locs]


{22.313, 33.4695, 25.1021, 12.5511, 6.56425}

obs = Delete[obs, locs]


{25, 36, 27, 8, 4}

Total[ (obs - exp)^2/exp]


3.31031

• How do you find the relevant internal functions? – Ian Jun 2 '15 at 14:06
• @Ian, you might want to see this. – J. M. will be back soon Jun 2 '15 at 14:09
• @Ian, i usually use the poor-man's one-liner , e.g, ?? **earson*, to search. If you are lucky you get a short list of results. If a function foo is ReadProtected, use ClearAttributes[foo, ReadProtected] then try ?? foo` to see the code. – kglr Jun 2 '15 at 14:30
• @kguler: I should preface my statements with "No documented" from now on - +1 for another enlightening internals answer. – ciao Jun 2 '15 at 22:07
• @ciao, thank you for the vote. – kglr Jun 3 '15 at 7:22