Setting tolerances for Monte Carlo tests: standard error, false rejections and bit-exact replay

A simulation test needs a justified acceptance range, with a confidence level, a standard error and any allowance for approximation error. A default tolerance such as 0.005 or 1%, hidden in a test helper, gives the caller little basis for judging whether it suits the sample size N. At small N it may reject a correct implementation too often; at large N it may miss a bias the test was meant to detect.

A same-seed rerun needs a separate check. If the tool promises bit-exact output for a fixed seed, version and platform, compare the outputs without a tolerance. Runs with different seeds can then be checked against the theoretical value using a statistical criterion. Mixing these checks can hide a changed bit pattern or make a statistical false rejection look like a broken RNG. The examples below use Python and numpy.

Sampling error and the acceptance range

A Monte Carlo estimate of RTP or hit rate is a random variable, so comparing it for equality with the theoretical value can fail because of sampling variation. One approach is to check whether the theoretical value falls within the estimated interval, or to write the comparison as an inequality:

abs(estimate - theory) <= z * SE + delta

SE is the standard error of the estimator. For a two-sided normal approximation, z comes from the confidence level: 1.96 for 95%. Delta is a separate allowance for approximation error, there to absorb the cases where the theoretical value itself is not exact: a paytable figure that was rounded, or a theoretical value that was itself produced by some other approximation. Its unit and its source belong next to the criterion, otherwise in six months nobody will know where it came from.

The NIST e-Handbook describes the relationship between a test and its corresponding confidence interval: in its normal-mean example with known population SD, reject the null hypothesis when its parameter value falls outside the interval, using the same z for the rejection threshold and interval half-width. The forms above can be compared in that way, but adding delta widens the acceptance range, so the original false rejection rate cannot simply be carried over. The handbook defines alpha as the risk of rejecting a true null hypothesis. That is the false rejection rate discussed here, and the observed rate needs to be assessed against the chosen criterion.

The SE of the sample mean is the population SD divided by sqrt(N), on the assumption that samples are independent and identically distributed with finite variance. For a Bernoulli estimator such as hit rate, SD is sqrt(p(1-p)), so with N = 10,000 and p = 0.25 the SE is about 0.0043. A few limits apply before any of this is used:

The sampling unit has to be genuinely independent. If several spins within one round are linked, as in free spins or cascading wins, you cannot treat every spin as an independent sample and take the single-spin SD over the square root of the total spin count as the SE; independent rounds may be suitable sampling units, but the statistic must still match the quantity being estimated.
Quadrupling N halves SE. That is the scaling of the standard error; the error on any particular run is not guaranteed to halve with it.
Rare-event hit rates do not always suit the normal approximation; with small p, Np needs to be large enough. When SD is estimated from a small sample, the assumptions of a t-based method need checking before using it. This post only deals with hit-rate estimates at large N and moderate p.

False rejection rates of a fixed tolerance at different N

Here the common pattern and the z times SE criterion are run side by side over 2,000 independent experiments, counting how often each one rejects. The true hit rate is 0.25, the fixed tolerance 0.005, and the seed is 20261101 throughout.

import numpy as np

P_TRUE = 0.25        # theoretical hit rate
FIXED_TOL = 0.005    # the common pattern: a default tolerance unrelated to N
Z = 1.96             # z for 95%
RUNS = 2_000         # independent experiments

rng = np.random.default_rng(20261101)

for n in (10_000, 1_000_000):
    se = np.sqrt(P_TRUE * (1 - P_TRUE) / n)
    fixed_rej = z_rej = 0
    for _ in range(RUNS):
        hits = rng.random(n) < P_TRUE
        err = abs(hits.mean() - P_TRUE)
        fixed_rej += err > FIXED_TOL
        z_rej += err > Z * se
    print(f"N={n:>9,}  SE={se:.5f}  fixed tol {FIXED_TOL} = {FIXED_TOL/se:5.2f} SE")
    print(f"   fixed-tol rejection {fixed_rej/RUNS:.3f}   z*SE rejection {z_rej/RUNS:.3f}   nominal alpha 0.050")

# real drift of 0.003: the fixed tolerance passes it at large N, z*SE catches it
n = 1_000_000
se = np.sqrt(P_TRUE * (1 - P_TRUE) / n)
hits = rng.random(n) < P_TRUE + 0.003
err = abs(hits.mean() - P_TRUE)
print(f"\nactual p=0.253, N={n:,}: error {err:.5f}, "
      f"fixed tol {'pass' if err <= FIXED_TOL else 'reject'}, "
      f"z*SE {'pass' if err <= Z * se else 'reject'}")

N=   10,000  SE=0.00433  fixed tol 0.005 =  1.15 SE
   fixed-tol rejection 0.255   z*SE rejection 0.047   nominal alpha 0.050
N=1,000,000  SE=0.00043  fixed tol 0.005 = 11.55 SE
   fixed-tol rejection 0.000   z*SE rejection 0.055   nominal alpha 0.050

actual p=0.253, N=1,000,000: error 0.00319, fixed tol pass, z*SE reject

The same 0.005 is worth 1.15 SE at N = 10,000, where the example rejects a correct implementation roughly once in four runs, close to P(|Z| > 1.15) of about 0.25 under the normal approximation. At N = 1,000,000 it is worth 11.55 SE: no correct runs in this batch are rejected, but a hit rate off by 0.003 also passes. The z times SE criterion detects that bias at this N and stays near the nominal 5% rejection rate at both sample sizes. These proportions are themselves estimates from 2,000 experiments: the same seed reproduces them exactly, and a different seed will shift them a little.

Making the tolerance a required argument guarantees that someone typed a number, and nothing more. The more useful habit is to record, beside each checkpoint, how z and SE were derived and where delta came from, so a reviewer has something to check against. Deterministic checkpoints also need separating out: comparing a closed-form expected value against a second formula involves a tolerance that comes from floating-point precision and has nothing to do with N, and it should not share a constant with the statistical ones.

Multiple checks accumulate false rejections

A 95% pass rate per check does not give a 95% pass rate per batch. With 20 independent checks at 5% false rejection each, the probability that at least one rejects is 1 - 0.95^20, about 64%. The simulation below uses the same hit-rate estimator, 20 checks of 10,000 samples per batch over 2,000 batches, then repeats with each check’s alpha divided by 20.

import numpy as np

K = 20            # 20 independent checks per batch
ALPHA = 0.05      # nominal false rejection rate per check
N = 10_000        # samples per check
P_TRUE = 0.25
BATCHES = 2_000   # repeat 2000 batches

rng = np.random.default_rng(20261101)
se = np.sqrt(P_TRUE * (1 - P_TRUE) / N)
z = 1.96          # for ALPHA = 0.05

batch_fail = 0
for _ in range(BATCHES):
    # 20 checks, each drawing N samples independently
    hits = rng.random((K, N)) < P_TRUE
    err = np.abs(hits.mean(axis=1) - P_TRUE)
    batch_fail += np.any(err > z * se)

print(f"theory 1 - (1-alpha)^K = 1 - 0.95^20 = {1 - (1 - ALPHA) ** K:.4f}")
print(f"simulated (batches with at least one rejection out of {BATCHES}) = {batch_fail / BATCHES:.4f}")

# switch each check to alpha/K (Bonferroni) and look at the batch rate again
z_adj = 3.023     # two-sided z for alpha/K = 0.0025
batch_fail = 0
for _ in range(BATCHES):
    hits = rng.random((K, N)) < P_TRUE
    err = np.abs(hits.mean(axis=1) - P_TRUE)
    batch_fail += np.any(err > z_adj * se)
print(f"with alpha/K per check, batch rejection = {batch_fail / BATCHES:.4f} (target <= {ALPHA})")

theory 1 - (1-alpha)^K = 1 - 0.95^20 = 0.6415
simulated (batches with at least one rejection out of 2000) = 0.6390
with alpha/K per check, batch rejection = 0.0520 (target <= 0.05)

Under those assumptions, a nightly CI batch with a fresh seed each night would fail on roughly two days in three over time. If rerunning becomes routine, a real bias may be dismissed along with the false alarms. Dividing each check’s alpha by the number of checks reduces the batch risk, provided each check has a valid false rejection rate. Normal approximations and finite simulations do not guarantee an observed proportion below the target. The wider thresholds also mean that detecting a bias of the same size requires more samples. Choose the acceptable bias and error risks before calculating N per check. Adjusting the tolerance until CI passes can obscure what the test is capable of detecting.

Bit-exact replay and statistical checks are different things

Bit-exact replay can be an acceptance contract: for the same seed, version and platform, require identical complete output without a tolerance. The example below compares the complete output with np.array_equal; the truncated hash is printed only for people to read, and matching hashes alone would not establish bitwise identity. The inputs, order of RNG draws and parallelisation must also be fixed. Changes to those conditions may change the output even when a statistical check still passes.

A statistical check changes the seed and compares the estimate with the theoretical value using an acceptance range. It does not require the hashes to match, nor does it require every seed to produce a different hash. Checking the mean and its interval tests a property of the mean; it does not validate the whole distribution. Each line of a paytable still needs checking on its own. A new seed usually changes the output, but nothing guarantees it, so the example compares only the estimate.

import hashlib
import numpy as np

P_TRUE, N, Z = 0.25, 100_000, 1.96
SE = np.sqrt(P_TRUE * (1 - P_TRUE) / N)

def simulate(seed):
    rng = np.random.default_rng(seed)
    return rng.random(N) < P_TRUE

def digest(hits):  # display only; equality uses np.array_equal
    return hashlib.sha256(np.packbits(hits).tobytes()).hexdigest()[:16]

# check one: same seed must be bit-identical, no tolerance
a, b = simulate(42), simulate(42)
print("same seed 42:", digest(a), digest(b), "identical" if np.array_equal(a, b) else "DIFFERENT")

# check two: new seed, bits usually differ, compare the estimate against the range
for seed in (42, 43, 44):
    hits = simulate(seed)
    err = abs(hits.mean() - P_TRUE)
    print(f"seed {seed}: hash {digest(hits)}  estimate {hits.mean():.5f}  "
          f"error {err:.5f}  {'pass' if err <= Z * SE else 'reject'} (threshold {Z * SE:.5f})")

same seed 42: 205382a815fa93c5 205382a815fa93c5 identical
seed 42: hash 205382a815fa93c5  estimate 0.24948  error 0.00052  pass (threshold 0.00268)
seed 43: hash 4b56bb7791cc4cf4  estimate 0.24974  error 0.00026  pass (threshold 0.00268)
seed 44: hash 72ba14e71e60cf15  estimate 0.25030  error 0.00030  pass (threshold 0.00268)

Mixing the two goes wrong in both directions. If the replay check is given a tolerance, a change to the draw order alters the hash while the estimate stays inside the range, the test stays green, and the mismatch surfaces when a replay is needed for certification. If a statistical check is treated as a replay check, a seed that lands in the 5% tail reads as a broken RNG, and someone goes off to fix the RNG. How seeds are fixed in tests, kept unpredictable in production and replayed for disputes is covered in the post on game RNG in production and tests; what certification actually examines is in Only 86 Million Ways to Shuffle: What Game RNG Certification Checks.

Self-checks must be able to fail

The last point is about CI rather than statistics. When an item in a precision self-check is expected to PASS and does not, it should be marked FAIL and the process should exit non-zero; a warning alone gives CI no failure status to act on. Alongside the normal cases, include one that is deliberately out of range, like the p = 0.253 case above, and confirm it really does report FAIL. Otherwise nobody knows whether that comparison is being executed at all.

Notes

Write the acceptance range as z times SE plus delta: z from the confidence level, SE varying with N, delta a separately documented approximation allowance.
SE = SD / sqrt(N) assumes independent, identically distributed samples; linked spins may call for independent rounds as sampling units, with a statistic defined for the intended estimate. This post covers hit rates at large N and moderate p only.
A fixed 0.005 tolerance rejected a correct implementation 25.5% of the time at N = 10,000 and passed a real 0.003 drift at N = 1,000,000; z times SE stayed near the nominal 5% at both sizes.
20 independent checks at alpha = 0.05 reject at least one in 64.15% of batches by theory and 63.9% in simulation; alpha / 20 per check brought the batch rate to 5.2%.
Under fixed execution conditions, the same seed should be checked against the complete output without a tolerance, and the truncated hashes in the example do not prove bitwise identity; a new seed is checked against z times SE without requiring matching hashes.
A self-check expected to PASS must exit non-zero when it fails, and a deliberately out-of-range case should confirm that it does.