Simulation asks a different question of a fit

A goodness-of-fit test asks whether a model describes the sample you have. A Monte Carlo model asks it to generate values you have never seen. These are not the same requirement, and a distribution can satisfy the first while failing the second.

Some simulation decisions depend on extreme quantiles, overflow probabilities or rare delays. The fit should be assessed against those output quantities. Different goodness-of-fit statistics weight discrepancies differently; none guarantees accurate extrapolation into a sparsely observed tail.

That tail is precisely the part of the distribution your sample knows least about.

How little the sample says about the tail

In a real fit of 500 observations (the full run), the two best-ranked models were an inverse Gaussian and a lognormal. Across the observed range they are nearly identical:

Quantile Empirical Inverse Gaussian Lognormal
0.50 3.3260 3.2595 3.2841
0.90 6.4504 6.6341 6.5473
0.95 8.0259 8.0596 7.9617
0.99 10.7939 11.4245 11.4908
0.999 16.4482 16.3888 17.3363

The fitted 95th percentiles are close to the empirical value. At the 99.9th percentile the two models differ by about 6%. With 500 observations, the empirical estimate interpolates between the two largest values and is highly uncertain.

Beyond the data the gap widens:

Threshold Empirical Inverse Gaussian Lognormal
P(X > 12) 0.0080 0.0076 0.0080
P(X > 15) 0.0040 0.0019 0.0024
P(X > 20) 0.0000 0.0002 0.0004

At 20 — above the sample maximum of 16.49 — one model says such an event is twice as likely as the other. Both passed every goodness-of-fit test. Neither is contradicted by the data, because the data are silent there.

A sample of 500 has only about five observations above its empirical 99th percentile. It has no observed exceedances of its maximum. Predictions beyond that point rely strongly on the fitted family and its assumptions; the absence of exceedances does not establish that their probability is zero.

What to do about it

Choose the tail deliberately. Since the data cannot settle it, the choice should be made on grounds you can state. Does the process have a physical maximum? Is there a mechanism that produces occasional extreme values — retries, escalations, compounding? A lognormal, a Weibull with shape below 1 and a gamma imply very different answers about rare events, and the difference is a modelling assumption, not a fitted result.

Check support before anything else. A model with unbounded support will eventually generate a value your process cannot produce. In a long run, “eventually” arrives. If a duration cannot exceed a shift length, either use a bounded family or truncate explicitly and document it.

Run several plausible candidates. Consider the information-criterion differences, diagnostics and support together, then run the simulation under more than one defensible model and compare the outputs that matter:

for name in ["inverse_gaussian", "lognormal", "inverse_gamma_3p"]:
    params = phi.sorted_distributions[name]["parameters"]
    # instantiate the distribution, draw the replications,
    # and record the output metric your decision depends on

If the answers agree, the choice between them was never important and you can stop worrying about it. If they diverge, that spread is your result. Reporting a single number from a single fitted distribution presents a modelling choice as a measurement.

Inspect tail discrepancies explicitly. AD and Q–Q plots can help, alongside uncertainty in the probabilities or quantiles of interest. No test is uniformly best for every alternative. See the test interpretation guide.

Consider not fitting at all. If you have enough observations and only need to reproduce the range you have seen, resampling the empirical distribution avoids every extrapolation problem above. Its limitation is exactly the mirror image: it can never generate a value larger than your maximum. Fit a parametric model when you need the tail beyond the data — and then be explicit that the tail is an assumption.

Fitting marginals is not fitting the process

Fitting each input separately and sampling independently removes dependence between inputs and across time. This can change aggregate variance and tail risk; the direction depends on the dependence structure and the output being modelled. Shared causes and serial patterns deserve explicit investigation.

No amount of care in fitting the individual distributions fixes this. Check for correlation and autocorrelation in the raw data before deciding that independent marginals are an adequate model.

A working protocol

  1. Fit, and read the top group rather than the top row.
  2. Discard anything whose support is wrong for the process.
  3. Compare the candidates’ behaviour at the quantiles your model consumes — not at the median.
  4. Run the simulation with each surviving candidate.
  5. Report the spread across them as part of the result.
  6. State separately what the tail assumption is and what evidence, if any, supports it.

The distribution is an input to the model, and an uncertain one. Treating the fit as settled because a test did not reject it moves that uncertainty out of sight rather than out of the answer.

Primary references