7
$\begingroup$

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[1];
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[[1]]
(* False *)
$\endgroup$
  • $\begingroup$ There is no mechanism I'm aware of to pull out the bins. $\endgroup$ – ciao Jun 2 '15 at 5:18
4
$\begingroup$

Statistics`GoodnessOfFitTestingDump`pearsonChiSquareTest

Key calculations for PearsonChiSquareTest are performed inside the function Statistics`GoodnessOfFitTestingDump`pearsonChiSquareTest.

enter image description here

SeedRandom[1];
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:

Statistics`GoodnessOfFitTestingDump`pearsonChiSquareTest[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:

Statistics`Library`DiscreteUnivariateDistributionQ[dist]

True

Bins

nbins = Statistics`GoodnessOfFitTestingDump`iOptimalBinCount[n]

13

bins = Quiet[N[Sort[Quantile[dist, N[Range[0, 1, 1/nbins], mp]], Less], mp]];
bins = bins + N[1/2, mp];
bins[[1]] = bins[[1]] - 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 = Statistics`GoodnessOfFitTestingDump`iCDF[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

$\endgroup$
  • $\begingroup$ How do you find the relevant internal functions? $\endgroup$ – Ian Jun 2 '15 at 14:06
  • 1
    $\begingroup$ @Ian, you might want to see this. $\endgroup$ – J. M. will be back soon Jun 2 '15 at 14:09
  • 1
    $\begingroup$ @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. $\endgroup$ – kglr Jun 2 '15 at 14:30
  • $\begingroup$ @kguler: I should preface my statements with "No documented" from now on - +1 for another enlightening internals answer. $\endgroup$ – ciao Jun 2 '15 at 22:07
  • $\begingroup$ @ciao, thank you for the vote. $\endgroup$ – kglr Jun 3 '15 at 7:22

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.