Inconsistent DistributionFitTest results?

Question

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.

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.