Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

In pondering this question, I used Mariana's function to generate the following data.

data = {356, 403, 49, 677, 109, 566, 111, 233, 189, 395, 72, 103, 394,
108, 255, 201, 197, 101, 112, 144, 262, 231, 171, 349, 522, 262, 
189, 128, 97, 188, 285, 459, 182, 220, 301, 154, 243, 250, 199, 
293, 141, 302, 64, 196, 106, 560, 115, 172, 54, 236, 183, 133, 218,
614, 111, 161, 310, 224, 134, 427, 130, 200, 380, 87, 430, 183, 
800, 368, 210, 221, 105, 104, 78, 213, 103, 586, 395, 312, 384, 
203, 141, 224, 107, 106, 172, 304, 141, 298, 250, 226, 268, 288, 
108, 116, 347, 123, 622, 135, 223, 229, 79, 74, 144, 88, 130, 284, 
272, 500, 310, 325, 247, 149, 612, 41, 100, 257, 229, 400, 486, 
142, 140, 136, 56, 411, 489, 83, 142, 59, 108, 264, 108, 160, 347, 
129, 137, 120, 100, 247, 117, 188, 121, 132, 316, 280, 336, 227, 
197, 156, 397, 144, 101, 317, 624, 171, 189, 72, 276, 261, 102, 92,
131, 384, 256, 87, 109, 390, 97, 62, 172, 311, 188, 506, 239, 269,
403, 356, 268, 397, 214, 202, 321, 148, 120, 169, 74, 75, 235, 
129, 90, 423, 514, 63, 233, 61, 82, 104, 167, 251, 198, 203, 316, 
309, 310, 305, 743, 334, 95, 169, 185, 1074, 126, 278, 343, 857, 
119, 80, 102, 92, 223, 151, 309, 127, 253, 346, 286, 240, 251, 413,
101, 158, 462, 77, 138, 333, 275, 223, 224, 123, 129, 251, 72, 
225, 174, 237, 530, 110, 295, 153, 136, 183, 137, 79, 182, 187, 
177, 152, 293, 165, 124, 118, 163, 154, 222, 111, 110, 67, 96, 269,
255, 190, 297, 72, 216, 129, 166, 83, 52, 252, 168, 82, 491, 208, 
427, 470, 462, 110, 365, 465, 135, 131, 165, 166, 420, 190, 511, 
928, 246, 349, 274, 184, 291, 145, 298, 470, 232, 302, 212, 182, 
209, 730, 106, 105, 761, 91, 124, 244, 351, 119, 462, 101, 262, 
233, 146, 512, 156, 138, 155, 76, 385, 168, 146, 430, 172, 208, 
121, 170, 271, 206, 120, 233, 210, 953, 353, 186, 199, 221, 272, 
494, 136, 292, 107, 265, 162, 235, 185, 214, 90, 167, 315, 238, 
109, 102, 425, 713, 149, 438, 41, 247, 233, 145, 268, 580, 174, 
115, 132, 99, 136, 140, 223, 149, 371, 520, 300, 301, 117, 69, 403,
143, 941, 107, 126, 234, 212, 139, 197, 558, 133, 45, 82, 91, 118,
554, 457, 340, 239, 600, 222, 136, 211, 182, 359, 171, 96, 161, 
68, 181, 118, 171, 226, 121, 309, 222, 149, 95, 304, 177, 204, 194,
123, 129, 126, 160, 353, 108, 249, 170, 326, 620, 83, 252, 104, 
134, 246, 154, 268, 152, 303, 143, 168, 422, 298, 186, 128, 97, 92,
316, 100, 182, 230, 198, 140, 217, 823, 371, 457, 122, 257, 207, 
53, 260, 112, 190, 66, 244, 267, 98, 210, 276, 189, 61, 107, 123, 
180, 93, 213, 207, 233, 155, 541, 339, 95, 314, 77, 314, 219, 609, 
354, 121, 208, 272, 244, 201, 134, 428, 45, 214, 254, 115, 223, 
145, 155, 287, 60, 138, 382, 132, 124, 218, 256, 255, 221, 142, 
246, 116, 184, 275, 147, 161, 378, 156, 149, 492, 91, 143, 325, 
181, 301, 275, 255, 240, 249, 578, 136, 177, 160, 107, 395, 151, 
233, 149, 386, 38, 214, 243, 188, 582, 513, 176, 234, 87, 70, 130, 
321, 123, 450, 125, 145, 594, 164, 600, 54, 335, 124, 310, 262, 
470, 442, 338, 219, 73, 951, 158, 229, 139, 129, 364, 257, 231, 
392, 468, 136, 157, 222, 108, 351, 306, 78, 121, 137, 347, 128, 
239, 219, 92, 259, 213, 98, 151, 170, 202, 446, 336, 293, 174, 183,
100, 345, 203, 194, 280, 330, 251, 335, 202, 198, 371, 399, 241, 
588, 527, 305, 621, 101, 124, 516, 311, 192, 228, 281, 127, 351, 
116, 468, 126, 155, 237, 282, 470, 427, 150, 80, 438, 232, 180, 
128, 482, 169, 224, 105, 362, 136, 135, 94, 137, 172, 292, 186, 91,
109, 144, 304, 184, 239, 285, 232, 89, 131, 376, 153, 298, 60, 97,
83, 583, 187, 338, 196, 75, 125, 161, 294, 115, 182, 51, 328, 232,
68, 339, 322, 171, 57, 331, 235, 113, 127, 176, 165, 240, 213, 
310, 96, 250, 171, 221, 140, 115, 145, 186, 343, 188, 146, 226, 
559, 103, 348, 272, 157, 156, 296, 218, 143, 306, 435, 150, 380, 
121, 163, 213, 283, 155, 290, 156, 372, 212, 172, 120, 336, 280, 
152, 101, 202, 325, 160, 98, 91, 259, 135, 209, 385, 210, 147, 214,
644, 102, 76, 576, 133, 52, 424, 187, 628, 421, 147, 211, 276, 
468, 592, 99, 391, 302, 191, 441, 164, 136, 223, 212, 101, 122, 
274, 198, 161, 648, 243, 210, 346, 330, 311, 123, 484, 183, 215, 
450, 255, 680, 532, 569, 102, 97, 151, 321, 151, 164, 198, 289, 
171, 103, 118, 172, 101, 340, 176, 206, 70, 233, 170, 190, 448, 
339, 387, 33, 239, 295, 200, 131, 322, 111, 516, 313, 365, 203, 85,
134, 134, 191, 228, 270, 125, 80, 145, 272, 229, 106, 151, 117, 
289, 120, 644, 140, 247, 133, 525, 232, 109, 243, 74, 152, 516, 
311, 179, 247, 191, 308, 355, 102, 598, 382, 153, 108, 77, 197, 
210, 200, 83, 86, 315, 304, 243, 329, 397, 282, 140, 578, 129, 211,
293, 219, 113, 471, 260, 160, 179, 341, 622, 311, 187, 175, 403, 
140, 239, 258, 193, 358, 83, 241, 320, 457, 111, 206, 96, 179, 152,
158, 574, 199, 217, 189, 663, 336, 388, 258, 351, 362, 369, 155, 
66, 230, 501, 247, 330, 383, 202, 567, 349, 117, 161, 524, 349, 
197, 162, 121, 1005, 343, 325, 255, 59, 303, 79, 203, 505, 337, 
607, 272, 170, 190, 129, 503, 780, 304, 243, 272, 146, 135, 689, 
105, 287, 406, 119, 58, 466, 90, 194, 111, 69, 113, 262, 145, 95, 
79, 93, 154, 272, 245, 238, 135, 65, 90, 209, 154, 455, 77};
plot1 = Histogram[data, Automatic, "PDF", Frame -> True, 
PlotLabel -> "Histogram PDF"]

Now I specify a candidate symbolic Gamma distribution and use DistributionFitTest

candidatedist = GammaDistribution[a, b, c, d];
fittestdatatable = 
DistributionFitTest[data, candidatedist, "TestDataTable", 
Method -> Automatic] 
DistributionFitTest[data, candidatedist, "TestConclusion"] 
dist2 = DistributionFitTest[data, candidatedist, "FittedDistribution"] 
pdf2 = PDF[dist2, x]
plot2 = Plot[pdf2, {x, 0, 1500}, PlotStyle -> Thick, Frame -> True, 
PlotLabel -> "dist2 PDF", MaxRecursion -> 2, PlotPoints -> 1000, 
PlotRange -> All];
plot3 = Show[plot1, plot2, 
PlotLabel -> "dist2 vs Histogram Comparison"]

The resulting output tells me that it failed to converge within 100 iterations and that the fit is rejected at the 5% level presumably due to the low P-value of 0.0409 via Pearson Chi Square. Yet the PDF/Histogram plot and the ProbabilityPlot both look good.

Q1. Is there a way to use something like MaxIterations to go beyond 100 ?

Q2. Did it truly "not converge" ?

Next, I use the very parameters found (a,b,c and d) by DistributionFitTest to define a specific numerical Gamma distribution and repeat the analysis.

a = 7.22667435871025;
b = 4.0245321518298836;
c = 0.5170201479394287;
d = 31.933431399237044;
dist3 = GammaDistribution[a, b, c, d]
DistributionFitTest[data, dist3, "TestDataTable", Method -> Automatic] 
DistributionFitTest[data, dist3, "TestConclusion"] 
pdf3 = PDF[dist3, x]
plot5 = Plot[pdf3, {x, 0, 1500}, PlotStyle -> Thick, Frame -> True, 
PlotLabel -> "dist3 PDF", MaxRecursion -> 2, PlotPoints -> 1000, 
PlotRange -> All];
plot6 = Show[plot1, plot5, 
PlotLabel -> "dist3 vs Histogram Comparison"]

Now Mathematica appears quite happy with the numerically defined Gamma distribution and produces a P-value of 0.9625 via Cramer-von Mises and reports that the distribution is not rejected at the 5% level. And as before (and as expected) the PDF/Histogram and ProbabilityPlots look good and identical to the previous plots.

Q3. Why is the same distribution first rejected at the 5% level and then subsequently not rejected at the 5% level? Which is the truth ?

Thanks for any insight.

share|improve this question

1 Answer 1

There is actually a fair bit going on here that can make this confusing. The critical thing is the difference between testing fit to a family of distributions compared to testing fit to a particular distribution. Let me demonstrate.

SeedRandom[23];
data = RandomVariate[NormalDistribution[1, 2], 100];

DistributionFitTest[data, NormalDistribution[mu, sigma]]

(* 0.088117 *)

DistributionFitTest[data, NormalDistribution[1, 2]]

(* 0.495841 *)

What is going on here? In the first case we are testing fit against all normal distributions (i.e. is the data normally distributed) and in the second case we are testing fit against a particular distribution (i.e. does the data follow a normal distribution with mean 1 and sd 2). It shouldn't come as a surprise that we get a better fit when testing against the distribution that the data actually came from.

Now things get much more interesting when we plug the "fitted distribution" back in to DistributionFitTest.

htd = DistributionFitTest[data, NormalDistribution[a, b], "HypothesisTestData"];

DistributionFitTest[data, htd["FittedDistribution"]]

(* 0.549637 *)

Whoa! What gives?

The general scheme for testing against a family of distributions is to test fit against the best fitting distribution in that family. But this naturally causes us to detect a better fit to the family than we should. Thus, when possible, a correction is applied to the p-value to adjust for the fact that we estimated the parameters prior to testing for fit.

When you inject the fitted distribution before hand DistributionFitTest doesn't know to correct the p-value and so the result is inflated.

Now for your particular example. As a rule DistributionFitTest attempts to fit the distribution family with the standard settings of FindDistributionParameters. This can fail and that is why you get the messages. For now, I'm not sure that there is a good workaround using DistributionFitTest itself. The difficulty lies in the fact that correction factors generally have to be hard coded for each distribution and parameter combination and often for each fitting method that can be used.

For what it's worth, most of what I'm saying here can be gleaned from the properties and relations examples for DistributionFitTest.

EDIT:

Incidentally the Lilliefors correction to the Kolmogorov-Smirnov test explains the motivation for all of this p-value correction business.

EDIT 2:

In order to better address the problem at hand I will borrow from my answer to this question. Since we can't use DistributionFitTest directly to do what we want we can always code up a simple chi-square test that corrects for the number of estimated parameters.

pearsonTest[obs_List, exp_List, cor_] /; Length[obs] == Length[exp] :=
  Block[{t}, t = Total[(obs - exp)^2/exp] // N;
  {Rule["chisqr", t], 
   Rule["p-val", 
    SurvivalFunction[ChiSquareDistribution[Length[exp] - 1 - cor], 
     t]]}]

n = Length[data];
nbins = Ceiling[2 n^(2/5)];
bins = Quantile[EstimatedDistribution[data, GammaDistribution[a, b, c, d]], 
          Range[0, 1, 1/nbins]] /. {\[Infinity] -> $MaxMachineNumber};

observed = HistogramList[data, {bins}][[2]];
expected = ConstantArray[n/nbins, nbins];

For the Gamma distribution you have given we are estimating 4 parameters from the data.

pearsonTest[observed, expected, 4]

(*{"chisqr" -> 41.024, "p-val" -> 0.0409676}*)

This suggests that the distribution may not fit particularly well. As can be seen from a quantile plot this is likely due to the divergence of the tails.

share|improve this answer
    
Thank you for your response but I can't reproduce/understand your results for 2 reasons. Firstly I don't reproduce your 0.257562 result. My version 9.0.1.0 returns 0.088117. Secondly, I get the following error message DistributionFitTest::rctnlndst: The argument htd[FittedDistribution] at position 2 should be a valid distribution or a rectangular array of real numbers with length greater than the dimension of the array. The dimensionality of the arguments at positions 1 and 2 must match. >> from your last line of code. –  Steve Mar 26 at 13:39
    
@Steve I get the same P value for the first answer as you do, and Andy's code is missing something like htd = DistributionFitTest[data, NormalDistribution[mu, sigma], "HypothesisTestData"] –  bobthechemist Mar 26 at 13:42
    
@Andy Ross, I can accept the fact that there are differences between families of distributions and a specific distribution but I think the larger question is the potential for a user to be misled. Regarding my original questions would you agree that the answers are 1)No, 2)No, 3)the truth is the Gamma distribution cannot be rejected at the 5% level ? Note that I'm not a Statistician or Mathematica power-user, just an Engineer needing to make decisions based on Mathematica's outputs. Thanks. –  Steve Mar 26 at 14:03
    
@Steve, sorry I was unsuccessfully trying to multitask. I've fixed the copy-paste errors. –  Andy Ross Mar 26 at 15:33
    
To test goodness of fit I recommend looking at a QuantilePlot against the fitted distribution you give. It appears that it fits pretty well until we get out into the right tail. –  Andy Ross Mar 26 at 15:39

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.