I have the following grouped data:

enter image description here

Can someone please tell me how to run a goodness of fit test to test if this data is from a normal distribution?

I would like to do this without estimating the parameters from the data.

I could only find PearsonChiSquaredTest, but it is for raw data, not grouped data.

  • $\begingroup$ This seems to be more a mathematical / statistics problem and not one for the Mathematica context. $\endgroup$
    – mgamer
    Commented Sep 9, 2019 at 15:00
  • 1
    $\begingroup$ Where did you find PearsonGoodnessOfFit? I see a PearsonChiSquareTest in version 12. You might be better off asking this question on CrossValidated and then come back to implement the process in Mathematica. Roughly you'll need to estimate the mean and standard deviation of grouped data and then construct a chisquare test accounting for the fact that you've had to estimate the mean and standard deviation. $\endgroup$
    – JimB
    Commented Sep 9, 2019 at 15:11
  • $\begingroup$ Sorry, I meant PearsoChiSquaredTest (I edited it now). My question was really whether there is a command in Mathematica to do this or not (so I did not intend it to be a theoretical question). This seems to be a natural question, so I was wondering if I missed a command or not. $\endgroup$ Commented Sep 9, 2019 at 15:27

1 Answer 1


You should still ask this on CrossValidated as I hope you'll get lectures about P-values. Here's one excellent lecture: Frank Harrell on P-values. (In short, one should not rely solely on P-values to make decisions.)

In the meantime here are two approaches to get approximate P-values. First define a function to get the chisquare value:

(* Some data *)
freq = {45, 55, 38, 27, 25, 10};
boundaries = {0, 10, 15, 20, 25, 30, 40};

(* Expand to allow for (nearly) the complete range of a normal distribution *)
freq = Join[{0}, freq, {0}]
(* {0, 45, 55, 38, 27, 25, 10, 0} *)
boundaries = Join[{-$MaxMachineNumber}, boundaries, {$MaxMachineNumber}]
(* {-1.79769*10^308, 0, 10, 15, 20, 25, 30, 40, 1.79769*10^308} *)

chisquare[freq_, boundaries_] := Module[{logL, midPoints, μ0, σ0, sol, d, expected, n},
  logL = Sum[freq[[i]]*
    Log[CDF[NormalDistribution[μ, σ], boundaries[[i + 1]]] - CDF[NormalDistribution[μ, σ], boundaries[[i]]]],
    {i, Length[freq]}];

  (* Get starting values for μ and σ *)
  midPoints = (boundaries[[2 ;; 7]] + boundaries[[3 ;; 8]])/2;
  (μ0 = midPoints.freq[[2 ;; 7]]/Total[freq]);
  (σ0 = Sqrt[((midPoints - μ0)^2).freq[[2 ;; 7]]/Total[freq]]);
  (* Find maximum likelihood estimates of μ and σ *)
  sol = FindMaximum[{logL, σ > 0}, {{μ, μ0}, {σ, σ0}}];

  n = Total[freq];
  d = NormalDistribution[μ, σ] /. sol[[2]];
  expected = Table[n (CDF[d, boundaries[[i + 1]]] - CDF[d, boundaries[[i]]]), {i, 1, 8}];
  {Total[(freq - expected)^2/expected], μ, σ} /. sol[[2]]]

t = chisquare[freq, boundaries]
(* {15.3559, 16.2655, 7.92905} *)

So the observed chisquare value is 15.3559. (The other two values are the respective maximum likelihood estimates of the mean and standard deviation from the grouped data.)

The approximate P-value using the chisquare distribution approximation is given by

1 - CDF[ChiSquareDistribution[Length[freq] - 3], t[[1]]]
(* 0.017662123393318407` *)

as there are 3 restrictions: (1) total count, (2) estimation of mean, and (3) estimation of standard deviation. (Note that I did not combine groups when expected counts are less than 5. One could and should add in that modification.)

A parametric bootstrap approach randomly samples from a normal distribution with mean and standard deviation equal to the maximum likelihood estimates. Then the chisquare statistic is calculated for each sample. The estimated P-value is proportion of bootstrap chisquare values greater than or equal to the observed chisquare value.

(* Now perform a parametric bootstrap to estimate P-value *)
d = NormalDistribution[t[[2]], t[[3]]];
n = Total[freq];
nBoot = 1000;
pValue = 0;
Do[x = RandomVariate[d, n];
 h = HistogramList[x, {boundaries}];
 c = chisquare[h[[2]], h[[1]]];
 If[c[[1]] >= t[[1]], pValue = pValue + 1],
 {i, nBoot}]
(pValue = pValue/nBoot) // N
(* 0.059` *)

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.