Statistics

core

Statistical analysis of simulation output.

Everything here works on plain sequences of numbers (one value per independent replication, or per Monte Carlo iteration). The confidence intervals assume independent, identically distributed observations; for autocorrelated output from a single long run use :func:batch_means first. Nothing here can guarantee the model itself is right - only that the sampling error of the reported estimates is quantified.

ReplicationAdvice dataclass

ReplicationAdvice(n: int, mean: float, half_width: float, relative_half_width: float, target_relative_half_width: float, required_n: int | None, sufficient: bool, note: str = '')

Answer to "have I run enough replications?" for one metric.

Uses the sequential approximation from Law (2015), Simulation Modeling and Analysis, sec. 9.4.1: the smallest n such that the t-based relative half-width, using the current variance estimate, is within target. It is an estimate: the variance itself is uncertain, so re-check after running the suggested number.

BatchMeans dataclass

BatchMeans(batch_means: FloatArray, summary: Summary, lag1_autocorrelation: float)

batches_look_independent property

batches_look_independent: bool

Heuristic: |lag-1 autocorrelation| within the ~95% band 2/sqrt(k) expected for k independent batch means. Passing does not prove independence.

Difference dataclass

Difference(mean_a: float, mean_b: float, difference: float, ci_low: float, ci_high: float, p_value: float, method: str, n_a: int, n_b: int, confidence: float)

Estimated difference b - a with a confidence interval.

significant property

significant: bool

True when the CI excludes zero (no correction for multiple comparisons).

ControlVariateEstimate dataclass

ControlVariateEstimate(mean: float, ci_low: float, ci_high: float, raw_mean: float, raw_ci_low: float, raw_ci_high: float, coefficients: dict[str, float], n: int, confidence: float)

Mean of y adjusted with control variates, and how much variance that removed.

variance_reduction property

variance_reduction: float

1 - (CV half-width / raw half-width)^2: share of estimator variance removed.

mean_ci

mean_ci(values: Sequence[float] | FloatArray, confidence: float = 0.95) -> tuple[float, float]

Student-t confidence interval for the mean (NaNs dropped).

bootstrap_ci

bootstrap_ci(values: Sequence[float] | FloatArray, statistic: Callable[..., Any] = np.mean, *, confidence: float = 0.95, n_resamples: int = 2000, seed: int = 0) -> tuple[float, float]

Percentile bootstrap interval for any statistic (e.g. median, p95).

statistic must accept a 2-D array and an axis keyword (NumPy reductions do), which lets resampling run vectorised.

proportion_ci

proportion_ci(successes: int, n: int, confidence: float = 0.95) -> tuple[float, float]

Wilson score interval for a probability estimated from successes / n.

convergence

convergence(values: Sequence[float] | FloatArray, confidence: float = 0.95) -> Convergence

Running mean and t half-width after each replication (for convergence plots).

batch_means

batch_means(observations: Sequence[float] | FloatArray, n_batches: int = 20, confidence: float = 0.95) -> BatchMeans

Non-overlapping batch means for a CI from one long (steady-state) run.

Observations in the same run are usually autocorrelated (one long wait follows another), so the naive t-interval is too narrow. Averaging over large batches reduces the correlation; check lag1_autocorrelation and use fewer, larger batches if it is high.

paired_difference

paired_difference(a: Sequence[float] | FloatArray, b: Sequence[float] | FloatArray, confidence: float = 0.95) -> Difference

Paired-t CI for b - a - appropriate with common random numbers.

welch_difference

welch_difference(a: Sequence[float] | FloatArray, b: Sequence[float] | FloatArray, confidence: float = 0.95) -> Difference

Welch (unequal-variance) t CI for mean(b) - mean(a) from independent samples.

adjust_p_values

adjust_p_values(p_values: Sequence[float], method: AdjustMethod = 'holm') -> list[float]

Adjust p-values for multiple comparisons; NaNs are kept and not counted.

• bonferroni - p * m; controls the family-wise error rate.
• holm - Holm's step-down; also controls the family-wise error rate but is uniformly more powerful than Bonferroni.
• bh - Benjamini-Hochberg; controls the false discovery rate (expected share of false positives among the "significant" results).

control_variate

control_variate(y: Sequence[float] | FloatArray, controls: Mapping[str, Sequence[float] | FloatArray], means: Mapping[str, float], *, confidence: float = 0.95) -> ControlVariateEstimate

Control-variate estimator of E[y] (Law 2015, sec. 11.4; regression form).

Each control x_k is an output whose expectation means[k] is known exactly (e.g. the average sampled service time). Regressing y on x - mean gives an unbiased intercept (the adjusted mean) with a standard error from the regression residuals, using n - 1 - q degrees of freedom for q controls. It helps only when the controls are correlated with y; with few replications the estimated coefficients cost some precision, so keep q small.