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.

The reason is where each one puts its attention. Goodness-of-fit statistics are dominated by the region containing most of the data — the centre. A simulation’s output is usually dominated by the opposite region: the worst percentile of a delay, the probability of exceeding a capacity, the size of the buffer that absorbs the bad days. Run a million replications and the answer that matters is decided by the few thousand draws from the far 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.5475
0.95 8.0259 8.0596 7.9621
0.99 10.7939 11.4245 11.4915
0.999 16.4482 16.3888 17.3377

Both are excellent through the 95th percentile. At the 99.9th percentile they differ by 6%, and there the empirical column is meaningless anyway — with n = 500 it is interpolating from a single observation.

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.

This is the honest position: with 500 observations you have roughly 5 above the 99th percentile, and no information at all beyond your maximum. Any statement your simulation makes about that region comes from the shape assumption of the fitted family, not from your sample.

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 the top candidates, not just the winner. This is the single most useful habit. Take every model within about 2 AIC of the leader, run the simulation with each, and compare the outputs you actually care about:

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.

Prefer Anderson–Darling among the tests. It weights the tails, which is the region your simulation consumes. A model selected on Kolmogorov–Smirnov alone has been selected on the part of the distribution that matters least here. The comparison of the three tests covers why.

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

A separate failure, and a more expensive one. Fitting each input independently and sampling them independently discards every relationship between them. If service time and arrival rate move together, if today’s delay predicts tomorrow’s, if two components fail from a shared cause — independent marginal draws will produce a simulation that is far too well behaved. Its aggregate variance will be understated, and the extreme scenarios that justify the model in the first place will be too rare.

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