Fitting is easier to trust once you have watched it work on data whose answer you already know. Generate 200 values from a gamma with shape 2, and see whether the ranking table picks the gamma out—and how close the runner-up gets. Add a location shift and see whether turning the shift box on recovers it.

Then change n and leave everything else alone. At n = 30 no goodness-of-fit test on this page can tell your gamma from a lognormal or a Weibull; every p-value is comfortable and every candidate survives. At n = 5000 the same tests separate them easily. Push further and something less comfortable happens, which the note beside the tests in section 4 takes up.

The recorded to setting rounds the values on the way out, which is what a stopwatch or a spreadsheet does to a real measurement. Generate at full precision and at one decimal place and compare the K–S statistic: rounding creates ties, and ties inflate both K–S and Anderson–Darling. The bootstrap p-values beside them hold up anyway: the resamples are drawn at full precision, without the rounding, and a calibration check shows that realistic recording granularity barely moves them. Rounding coarse enough to matter tends to break the fit itself first—rounding small positive values to whole numbers creates zeros that no gamma, lognormal, or Weibull can absorb.

The two Johnson entries sit under their own heading in the menu because they are the only ones that cannot also be fitted. That is what they are for. Every other choice here rewards the ranking table for recovering a distribution that is genuinely on its menu, which is a flattering test; real data is never on anybody's menu. A Johnson sample gives the table nothing to recover; what you see is how well the best available approximation does—which is the situation you are actually in every time you fit a distribution to something you measured.

The Johnson system reaches almost any combination of skewness and kurtosis by putting a standard normal Z through a monotone transformation. The two entries here are different distributions, not one distribution with different starting numbers: they use different transformations, and no choice of parameters turns either into the other.

• Sᵇ is bounded: a logistic transformation, X = ξ + λ / (1 + exp(−(Z−γ)/δ)), which lives on the finite interval from ξ to ξ + λ and cannot leave it whatever γ and δ are. Reach for it when the quantity has real limits (a proportion, a fill level, a task with a floor and a ceiling).
• Sᵘ is unbounded: a hyperbolic-sine transformation, X = ξ + λ sinh((Z−γ)/δ), whose support is the whole real line. It is the member with genuinely heavy tails, which is why a large enough sample from it contains a negative value or two, and why the positive-only candidates then refuse until you turn the shift on.

Johnson's system has two more members, and between them the four cover the plane: the lognormal is Sᴱ, and the normal itself is Sₙ. Which one a dataset needs is settled by where its skewness and kurtosis fall, with the lognormal line acting as the border between the Sᵇ and Sᵘ regions. That is why fitting a Johnson distribution is a different kind of job from fitting anything else on this page, and part of why it is not offered here.

With new seed each time ticked, pressing Generate again draws a different sample from the same distribution. That is the quickest way to see how much of what the ranking table tells you is the distribution and how much is the luck of one sample: at n = 50, press it five times and watch the winner change.

Log-likelihood asks how probable your data is under each fitted model. It always improves when a distribution has more parameters to spend, and so on its own it would pick the most flexible candidate every time, whether or not the extra flexibility described anything real.

AIC and BIC are the log-likelihood with a price attached to each parameter: AIC = 2k − 2 lnL and BIC = k ln n − 2 lnL, both smaller-is-better. Because ln n is larger than 2 for any sample past seven observations, BIC charges more per parameter than AIC and so leans toward simpler models, and it leans harder the more data you have. When they disagree, the disagreement is telling you the extra parameter is doing something worth about one or two log-likelihood units—real, but marginal.

Square error is the criterion Arena's Input Analyzer ranks by. It lays the fitted density over the histogram and adds up the squared gaps between the observed relative frequency in each bin and the probability the fitted model puts in that bin. Two things are wrong with using it to choose:

• It depends on the binning. The same data and the same two candidates can change places when you drag the bin slider because the criterion is computed on the histogram rather than on the data. Try it: sort by square error and then move the slider in section 4.
• It does not charge for parameters. A three-parameter distribution can always track a histogram at least as closely as a one-parameter distribution can, and so square error rewards flexibility for its own sake—the same flaw as raw log-likelihood, without the compensating property of being a likelihood.

So use AIC or BIC to choose, use the Q–Q plot to check that the choice is sensible, and read square error only to understand what Arena would have told you.

Maximum likelihood picks the parameter values under which the data you actually collected would have been the most probable. For a few distributions that works out to a formula you can do by hand—the exponential's fitted mean is just the sample mean, and the normal's are the sample mean and standard deviation. For the rest there is no such formula, and the estimate is found by numerical search: the gamma and beta fits solve their likelihood equations by Newton's method, the Weibull's shape by bisection, and the triangular by Nelder–Mead, a search that needs no derivatives at all.

A location shift γ slides a distribution sideways so it starts somewhere other than zero. It applies only to the candidates that begin at zero—exponential, gamma, Erlang, Weibull, lognormal—because the others already carry a location parameter of their own and a second one would be redundant. Turn the shift on when the process has an unavoidable minimum: a fixed setup time, a mandatory inspection, a walk down a corridor that cannot be done in zero seconds. Turn it off otherwise because it costs a parameter and it is estimated from the smallest value you happened to see.

The chi-square test compares the count that landed in each histogram bin with the count the fitted distribution predicts for that same bin: χ² = Σ (O − E)² / E, referred to a chi-square distribution with (bins − 1 − parameters estimated) degrees of freedom. Every ingredient depends on where you put the bin edges, and so the statistic and its p-value are properties of your binning as much as of your data. Drag the slider and watch: the same data and the same fitted model will hand you a comfortable p-value at one bin count and an uncomfortable one a few notches away.

That is not a bug in the test, and it is not licence to shop for the bin count that gives the answer you want. It is the reason chi-square is the weakest of the three tests here, and the reason the other two do not bin at all. Too few bins and the test cannot see a misfit inside a wide bin; too many and the expected counts get so small that the chi-square approximation stops holding.

That second failure is handled rather than merely flagged. Before the statistic is computed, adjacent bins are pooled left to right until every group expects at least 5 observations, which is the long-standing rule of thumb for the approximation to be trustworthy. The degrees of freedom then count the pooled groups, not the bins you see, and the table below shows which bins ended up together. On a right-skewed fit this almost always merges the thin far-right bins into one group, which is exactly the region a histogram cannot support on its own. (No continuity correction is applied. Yates' correction belongs to 2×2 tables and would be out of place here.)

So the slider moves the p-value for two reasons at once: it changes which observations share a bin, and it changes how many groups survive pooling. Push the bin count high enough and the group count stops rising because there is only so much data to spread.

Chi-square does correct for having fitted the distribution to the same data, in the only place it can: it subtracts one degree of freedom per estimated parameter. Two details of that subtraction are worth knowing, and the readout applies both.

Not every estimated quantity costs a degree of freedom. A uniform's endpoints are the sample minimum and maximum; a beta's interval is read off the same two numbers; a shifted exponential's threshold is the sample minimum. An extreme order statistic converges to what it is estimating at rate 1/n instead of the usual 1/√n, and so asymptotically it behaves like a constant that was known all along and absorbs none of the discrepancy the correction exists to account for. Subtracting for it would throw away power for nothing; the degrees of freedom therefore count only the regular parameters. AIC and BIC keep the full count because their penalty is a charge for complexity rather than a distributional correction, and a number read off the data is a number the data paid for however fast it converged.

And the correction is exact only when the parameters were estimated from the binned counts, which is not what happens here—every fit on this page uses the ungrouped data. Chernoff and Lehmann showed in 1954 that the true distribution then sits somewhere between the two chi-squares you would get by subtracting nothing and by subtracting everything; the conventional df is an approximation whichever way you count. That is why the chi-square statistic is bootstrapped here too. The bootstrap re-bins and re-fits each resample exactly as your own sample was handled, and so its p-value needs no convention at all—and comparing it with the table column tells you how much the convention cost you on your data.

A bootstrap answers a question about a statistic when nobody can write down its distribution: it manufactures that distribution by simulation instead of deriving it on paper.

The question here is: if the fitted distribution really were the truth, and someone fitted it to a fresh sample exactly the way you fitted yours, how large would D usually come out? There is no formula because the answer depends on which distribution was fitted and on how many parameters were estimated. So the page simulates it:

• draw n fresh values from the fitted distribution, with n the same as your sample size;
• refit the parameters to those new values, by the same method;
• compute D for the new sample against its own new fit;
• repeat a couple of thousand times.

Those couple of thousand values of D are the null distribution: the range D lands in when the model is right. Your D is then read against them, and the p-value is just the fraction that came out at least as large. The refitting step is what makes it work—a version that reused your fitted parameters instead of estimating fresh ones would reproduce the table's error exactly.

It is not subsampling: your sample is never split, thinned, or reduced, and n is the same in every replication as in your data. Nor is it a fix for the tests being too powerful at large n—it does not touch your data or your statistic at all, only the yardstick the statistic is measured against. (The word parametric in "parametric bootstrap" marks the difference from the more familiar kind, which resamples your own observations with replacement. Here the resamples come from the fitted distribution, which is what makes them samples from the null hypothesis rather than from your data.)

The published Kolmogorov–Smirnov table—the one you can use by hand, and the one most software reports—answers a different question from the one you are asking. It gives the distribution of D when the distribution being tested against was fully specified before you saw the data: the family named, the parameters written down, everything. That distribution of D has the convenient property of not depending on which distribution you named, which is why one short table can serve every test.

You did not do that. You estimated the parameters from the very sample you are now testing, which tilts the fitted curve toward your data. The curve passes closer to the empirical CDF than a curve chosen in advance would have, so D comes out smaller than the table expects. The consequence runs in the direction people usually guess wrong: using the table anyway makes the test too forgiving, not too strict. It accepts fits it should reject. The table p-value shown beside the bootstrap one is almost always the larger of the two, and the gap between them is the size of the error.

There are published tables for a few special cases—Lilliefors' for a fitted normal and a fitted exponential, Stephens' for several more—but each is specific to one distribution and one way of estimating it. There is no such table for a fitted Erlang, or a fitted triangular, or a fitted gamma with a location shift. The bootstrap builds the right one for whatever you actually did.

Anderson–Darling is the same idea as K–S with the differences weighted to count more heavily in the tails, and it carries exactly the same caveat. For a queueing model the tails are usually what matter because the long service times are what build the queue; prefer A–D to K–S when the two disagree.

A goodness-of-fit test asks whether your data could plausibly have come from the fitted distribution. Once you have enough data, the honest answer is almost always no. Real service times are not exactly gamma: they are rounded to the nearest tenth of a minute, they drift between the morning and afternoon operator, they have a hard floor at the setup time, and one customer in fifty asks a question. None of that is visible at n = 30. All of it is visible at n = 5000, and every candidate gets rejected.

That is the test working, not failing—but it does mean a p-value stops being the useful question once n is large. What matters then is whether the departure is big enough to change the decision the simulation is for, and a p-value cannot tell you that. The Q–Q plot can: a fit that misses the upper tail by 10% will build the wrong queue, and a fit that is "wrong" only in the third decimal of the body will not.

The opposite trap sits at the other end. At n = 30 nothing gets rejected, and that is not evidence that your distribution is right—it is evidence that 30 observations cannot tell a gamma from a lognormal. "Not rejected" never means "correct"; at small n it usually means "not enough data to say".

Generate a sample from a known distribution in section 1 and step n from 30 to 5000, changing nothing else. You will watch both halves of this happen to a fit that was right the whole time.

Sort your n observations. The i-th smallest is an estimate of the quantile at probability (i − 0.5)/n; the plot puts it on the vertical axis against the fitted distribution's quantile at that same probability on the horizontal axis. A perfect fit puts every point on the 45° line.

The shape of the departure tells you what is wrong. Points bending above the line at the right-hand end mean your data has a heavier upper tail than the fitted distribution—the big values are bigger than the model expects, which for a queueing model is usually the departure that matters. An S-shape means the fitted distribution has the wrong skewness. A few stray points at the very ends, with everything else on the line, is normal because the extreme order statistics are the noisiest. Do not reject a fit over the last point alone.

This is the plot to fall back on when a p-value stops being useful, which at a large n it will.

The P–P plot compares probabilities where the Q–Q plot compares values: it puts the fitted F(x₀₋ₖ₋) against (i − 0.5)/n. Because probabilities are squeezed into [0,1], the tails get compressed and the middle gets stretched—exactly the opposite of the Q–Q plot's emphasis, where a handful of extreme values occupy most of the axis.

So use the P–P plot to judge the body of the distribution, where most of your observations are and where most of the probability lives, and the Q–Q plot to judge the tails. A fit that looks perfect on one and poor on the other is telling you where its problem is, which is more useful than either plot alone.

Every tool below describes the same fitted distribution, but they do not agree on what to call its parameters, and the disagreements are the kind that produce a model that runs fine and is quietly wrong:

• Exponential. Arena's EXPO and Simio's Random.Exponential take the mean. AnyLogic's exponential() and R's rexp() take the rate, which is one over the mean. Handing a rate to a function that wants a mean is the single most common input-modeling error, and nothing will warn you.
• Lognormal. Arena's LOGN takes the mean and standard deviation of the lognormal variable itself. R, MATLAB, and scipy take μ and σ of the underlying normal, that is, of the logarithm. The conversion is shown with the parameters in section 3, and the expressions below already have it applied.
• Weibull and gamma. Arena writes the scale first and the shape second, WEIB(scale, shape); R and scipy take the shape first. Both orders appear below, each correct for its own tool.
• The location shift. Arena has no shift built into its distributions, and so a shifted fit is written as an added constant, as in 0.5 + EXPO(2.3). Simio and scipy take it as an explicit offset argument.

The R and Python lines reproduce the fit and the goodness-of-fit test from your raw data rather than hard-coding the fitted numbers, and so running them is a genuine independent check of what this page told you—not a restatement of it.