api.sim

sim_runlength_plan(means, ciHalfWidth, samplesUsed, varargin)

PLAN = SIM_RUNLENGTH_PLAN(MEANS, CIHALFWIDTH, SAMPLESUSED)

How long a simulation run should have been, from the one it already did.

A batch-means half-width H at confidence 1-alpha over a run of N samples pins the ASYMPTOTIC variance of the estimator,

sigma^2 = (H/z)^2 N, z = Phi^-1((1+confidence)/2),

and that is the quantity a run length is planned from – NOT the stationary variance, which on M/M/1 differs from it by a factor blowing up like (1-rho)^-2. SIM_RUNLENGTH then turns it into the sample count that reaches a requested RELATIVE precision.

MEANS and CIHALFWIDTH are matrices of the same shape, one entry per (station, class); an entry with a non-positive mean or half-width is left NaN, since there is nothing to plan from there.

Options: ‘relprecision’ (default 0.05), ‘confidence’ (default 0.95, and it must be the level the half-widths were computed at).

Returns a struct with fields relprecision, confidence, samplesUsed, asymptoticVariance and requiredSamples.

Reference: W. Whitt (1989). Planning queueing simulations. Management Science 35(11), 1341-1366.

See also SIM_RUNLENGTH, SIM_ASYMVAR_CTMC, SIM_ASYMVAR_MM1.

sim_runlength(mean_, asymVar, varargin)

SIM_RUNLENGTH Run length for a steady-state estimate of a given precision.

RESULT = SIM_RUNLENGTH(MEAN, ASYMVAR) returns the simulated time needed to estimate a steady-state mean to within 5% at 95% confidence, given the ASYMPTOTIC VARIANCE of the estimator.

A time average over [0,t] has standard error sqrt(sigma^2/t), so a two-sided interval of half-width z sqrt(sigma^2/t) reaches relative precision eps when

t* = (z/eps)^2 sigma^2 / mean^2.

THE POINT OF THE FORMULA is that everything expensive sits in sigma^2/mean^2, the squared coefficient of variation of the TIME AVERAGE rather than of the process. Halving the tolerance quadruples the run.

Options: ‘relPrecision’ (default 0.05), ‘confidence’ (default 0.95), ‘runLength’ (an actual run length, to report the precision it buys).

Returns a struct with fields requiredRunLength, z and, when runLength is given, halfWidth and achievedRelPrecision.

Reference: W. Whitt (1989). Planning queueing simulations. Management Science 35(11), 1341-1366.

See also SIM_ASYMVAR_MM1, SIM_ASYMVAR_CTMC.

sim_asymvar_mm1(lambda, mu)

SIM_ASYMVAR_MM1 Asymptotic variance of the M/M/1 number-in-system process.

THE QUANTITY THAT MATTERS FOR PLANNING is not the variance of the process but its ASYMPTOTIC VARIANCE, sigma^2 = lim t Var(time-average over [0,t]), which is twice the integral of the autocovariance. It is what says how long a run must be, because the time average of a positively correlated process converges at rate sigma^2/t, not at Var(X)/t.

For M/M/1 with utilization rho,: E[N] = rho/(1-rho), Var(N) = rho/(1-rho)^2, sigma^2 = 2 rho (1+rho)/(mu (1-rho)^4).

The FOURTH power is the whole story: the variance of the process grows like (1-rho)^-2, but the run length needed to average it away grows like (1-rho)^-4 over the squared mean, i.e. like (1-rho)^-2. A queue at rho = 0.9 needs about 100 times the run of one at rho = 0.

Returns a struct with fields mean, variance, asymptoticVariance and relaxationTime (sigma^2/Var, the correlation time scale).

Reference: W. Whitt (1989). Planning queueing simulations. Management Science 35(11), 1341-1366; J. Abate, W. Whitt (1988). The correlation functions of RBM and M/M/1. Stochastic Models 4(2), 315-359.

See also SIM_RUNLENGTH, SIM_ASYMVAR_CTMC.

sim_asymvar_ctmc(A, f, pi)

SIM_ASYMVAR_CTMC Asymptotic variance of a reward on a CTMC.

RESULT = SIM_ASYMVAR_CTMC(A, F) returns sigma^2 = 2 sum_x pi(x) g(x) d(x), where g = F - E_pi[F] and d solves A d = -g: the DEVIATION vector, the accumulated future excess reward started from each state.

This is the general form of the quantity SIM_ASYMVAR_MM1 gives in closed form for M/M/1, and it is what run-length planning needs for any model LINE can build a generator for. A alone is singular: a constant may be added to d without changing the variance, so one normalization is needed to pin it.

Returns a struct with fields mean, variance, asymptoticVariance and deviation.

Reference: W. Whitt (1989). Planning queueing simulations. Management Science 35(11), 1341-1366.

See also SIM_RUNLENGTH, SIM_ASYMVAR_MM1, CTMC_SOLVE.

sim_vonneumann(x, alpha)

OA_VONNEUMANN Von Neumann ratio test for randomness of a sequence.

RESULT = OA_VONNEUMANN(X) applies the two-sided randomness test of von Neumann (1941) to the sequence X, at the default significance level 0.05.

RESULT = OA_VONNEUMANN(X, ALPHA) uses significance level ALPHA.

The statistic is the ratio of the mean square successive difference to the variance,

ratio = sum_{i=1}^{b-1} (x_{i+1}-x_i)^2 / sum_{i=1}^{b} (x_i - xbar)^2,

with b = numel(X). Under the null hypothesis that X is i.i.d. normal the ratio has mean 2 and variance 4(b-2)/((b-1)(b+1)), and (ratio-2)/sd is asymptotically standard normal, so the two-sided p-value is 2(1 - Phi(|z|)). Serial correlation of either sign moves the ratio away from 2: positive correlation shrinks the successive differences and pushes the ratio below 2, negative correlation pushes it above.

The null mean and variance above were confirmed by Monte Carlo over b = 10, 16, 24, 32, 50 to within 0.3%.

Returns a struct with fields:: ratio - The von Neumann ratio zscore - Standardized statistic (ratio-2)/sd pvalue - Two-sided p-value reject - true when pvalue < ALPHA, i.e. randomness is rejected nobs - Number of observations b analyzer - Identifier string

Examples

sim_vonneumann(randn(1,50)).reject % false, i.i.d. input sim_vonneumann(cumsum(randn(1,50))).reject % true, random walk

Reference: J. von Neumann, “Distribution of the Ratio of the Mean Square Successive Difference to the Variance”, Ann. Math. Statist. 12(4), 1941; L. C. Young, “Randomness in Ordered Sequences”, Ann. Math. Statist. 12, 1941.

See also OA_SHAPIROWILK, OA_FQUEST, OA_FIRQUEST

sim_tinv(p, nu)

OA_TINV Quantile function of Student’s t distribution.

T = OA_TINV(P, NU) returns the P-quantile of the t distribution with NU degrees of freedom. P is a scalar or array, NU a positive scalar.

Implemented on the incomplete beta inverse, which is in base MATLAB, so the output-analysis routines in this folder do not pull in the Statistics and Machine Learning Toolbox. The identity used is

P(|T| > t) = betainc(nu/(nu+t^2), nu/2, 1/2),

inverted for the two-sided tail 2(1-P) and mapped back with: t = sqrt(nu (1-z)/z), z = betaincinv(2(1-P), nu/2, 1/2).

Agreement with the toolbox TINV is to within 1e-10 over the range used here.

Examples

sim_tinv(0.975, 10) % 2.2281 sim_tinv(0.975, 29) % 2.0452

See also OA_NORMINV, OA_TINV

sim_fquest(Y, p, alpha, options)

OA_FQUEST Fixed-sample-size confidence interval for a steady-state quantile.

RESULT = OA_FQUEST(Y, P) returns a point estimate and a 95% confidence interval for the P-quantile of the steady-state marginal distribution of the simulation output process whose single sample path is Y. Y has arbitrary fixed length; no sequential control of the run length is needed.

RESULT = OA_FQUEST(Y, P, ALPHA) uses nominal coverage 1-ALPHA.

RESULT = OA_FQUEST(Y, P, ALPHA, OPTIONS) overrides the procedure constants. OPTIONS is a struct whose recognized fields are

b0 - initial batch count for the warmup stage (default 50) m0 - initial batch size for the warmup stage (default 500) s - descending batch counts for the test stages (default [32 24 16 10]) beta - significance level of the stage tests (default 0.30) eta - decay coefficient of the warmup significance (default 0.2) theta - decay exponent of the warmup significance (default 2.3) weight - constant STS weight function (default sqrt(12)) force - deliver a heuristic interval when a test fails (default true)

The procedure is FQUEST. It has four blocks:

Warmup. Starting from B = b0 and M = m0 it computes the B signed STS areas of the batched quantile process and tests them for randomness with von Neumann’s ratio at the decaying significance beta*exp(-eta*(l-1)^theta) on iteration l, growing M by sqrt(2) whenever the test rejects. Passing the test means the areas are approximately independent, so any initialization bias is confined to the first batch.

Truncation. The first batch is deleted, which is the entire warmup treatment; there is no separate transient detector.

Batch-count selection. With B stepping down through s and M = floor(N*/B), four tests must pass in order: von Neumann and Shapiro-Wilk on the signed areas, then von Neumann and Shapiro-Wilk on the batched quantile estimators. These check the asymptotic properties the interval rests on, namely that both the areas and the batch quantiles behave like independent normal variates. B only ever decreases, and a failure at B = 10 ends the stage.

Delivery. When all four tests pass the interval is
ytilde_p(n*) +- t_{1-alpha/2, 2b-1} sqrt(V_p(w;b,m)/n*)

with V_p the combined variance-parameter estimator of OA_STS_QUANTILE_AREAS. Otherwise the sample was too small, RESULT.heuristic is true, and with OPTIONS.force the interval returned is the union of the wider of the two single-component intervals and Willink’s skewness- and correlation-adjusted asymmetric interval, both built from the batched quantile estimators.

Returns a struct with fields:

estimate - Full-sample empirical P-quantile of the truncated path lower,upper - Confidence interval endpoints halfwidth - Half-width, (upper-lower)/2, which the asymmetric heuristic

interval attains only on average

b, m, n - Final batch count, batch size and number of observations used truncated - Number of observations deleted from the front of Y Ap, Np, Vp - The three variance-parameter estimators at the final b and m heuristic - true when a test failed and the interval is not asymptotically

justified

warnings - Cell array of diagnostic strings, empty on a clean run analyzer - Identifier string

Coverage was measured on the article’s own test bed, the waiting-time process of an M/M/1 queue with lambda = 0.8, mu = 1 started with 113 jobs in system, over 500 independent replications at N = 200000, giving a standard error near 1%: 95.2% at p = 0.5, 96.2% at p = 0.9 and 95.6% at p = 0.99 against a nominal 95%. The delivered half-width exceeds the empirically needed one by factors of 1.16, 1.24 and 1.99 respectively, so the interval is conservative and increasingly so into the tail, consistent with the half-widths the article reports. On i.i.d. Exp(1) data, where sigma_p^2 = p(1-p)/f(y_p)^2 is exact, both A_p and N_p are unbiased to within 7%.

Two properties are worth knowing. Coverage is not monotone in the sample size over this range: a larger sample passes the stage tests more often and so reaches the conservative fallback less. And a substantial fraction of runs takes that fallback at all, 22% to 65% here and rising with p, so a delivered interval may well be the heuristic one; RESULT.heuristic says which. Fewer than about 100 replications cannot resolve a two-point difference in coverage, so do not read a small experiment as a defect.

Applicability is a condition on the output process, not on the model that produced it. The theory needs geometric moment contraction (Wu 2005), which holds for ARMA series, a broad class of short-range-dependent linear and nonlinear processes, many Markov chains, and was proved for M/M/1 and non-heavy-tailed G/G/1 waiting times by Dingec et al. (2022); a density that is positive and differentiable at the quantile of interest; short-range dependence and an FCLT for the indicator process. M/M/1 is only the validation bed, chosen because its exact quantiles are known.

Two practical exclusions follow. Do NOT use this on integer-valued output such as a queue length: the marginal has no density, the density-regularity condition fails, and the batched quantile has no Bahadur representation. Use it on continuous output, that is response, waiting and sojourn times. And heavy-tailed service, which can break geometric moment contraction and induce long-range dependence, is outside the theory.

Examples

Y = exprnd(1, 200000, 1); r = sim_fquest(Y, 0.9, 0.05); [r.lower r.estimate r.upper] % brackets -log(0.1) = 2.3026

Reference: A. Lolos, C. Alexopoulos, D. Goldsman, K. D. Dingec, A. C. Mokashi, J. R. Wilson, “A Fixed-Sample-Size Method for Estimating Steady-State Quantiles”, Proc. Winter Simulation Conference, 2023.

See also OA_FIRQUEST, OA_STS_QUANTILE_AREAS, OA_VONNEUMANN, OA_SHAPIROWILK

sim_firquest(Y, p, alpha, options)

OA_FIRQUEST Fixed-sample-size quantile interval from independent replications.

RESULT = OA_FIRQUEST(Y, P) returns a point estimate and a 95% confidence interval for the P-quantile of the steady-state marginal distribution, from R independent replications of equal length. Y is an N-by-R matrix whose columns are the replicate sample paths, or a cell array of R equal-length vectors.

RESULT = OA_FIRQUEST(Y, P, ALPHA) uses nominal coverage 1-ALPHA.

RESULT = OA_FIRQUEST(Y, P, ALPHA, OPTIONS) overrides the procedure constants; see OA_FQUEST for the recognized fields. Two defaults differ from FQUEST: b0 = 25 rather than 50, and s is chosen from R rather than fixed, because the stage tests act on the R*B pooled statistics and so need fewer batches per replication. Supplying s explicitly overrides that choice.

The procedure is FIRQUEST, the replicated counterpart of FQUEST. It differs from OA_FQUEST in four places:

The warmup randomness test runs independently on each replicate path, and the batch size it settles on may differ between replications.

Truncation removes the largest of those batch sizes from the front of every replication, not just from one path. This is more aggressive than FQUEST on purpose: an untruncated transient common to all replications biases every replicate estimate the same way, and averaging cannot remove it.

The four stage tests act on the R*B signed areas and R*B replicate batched quantile estimators pooled across replications, with the same B in every replication and at least one batch from each.

The delivered interval is
ytilde_p(N*) +- t_{1-alpha/2, 2Rb-1} sqrt(Vtilde_p(w;R,b,m)/N*),

N* = R*B*M, with the pooled combined variance-parameter estimator
A_p(w;R,b,m) = (Rb)^{-1} sum_j A_p(w;j,m)^2 Ntilde_p(R,b,m) = m (Rb-1)^{-1} sum_j (yhat_p(j,m) - ytilde_p(N*))^2 Vtilde_p = [Rb A_p + (Rb-1) Ntilde_p] / (2Rb-1).

The heuristic fallback drops FQUEST’s residual-autocorrelation correction, since the pooled batch quantiles come from independent paths.

The default batch counts s as a function of R are the article’s:: R = 2 -> [14 11 8 5], R = 3 -> [10 8 6 4], R = 4 -> [6 5 4 3], 5..9 -> [5 4 3 2], 10..16 -> [4 3 2 1], 17..22 -> [3 2 1], 23..32 -> [2 1], R >= 33 -> [1].
Returns a struct with the fields OA_FQUEST returns, plus: R - Number of replications truncated - Observations deleted from the front of every replication

and with b and m the per-replication batch count and batch size, so that n = R*b*m.

Independent replications shorten the correlation the estimator has to fight, and they parallelize, but they reintroduce initialization bias in every path, so a short run length per replication is worse here than in OA_FQUEST. The article reports slight undercoverage at P = 0.99 when the total sample is under 500000, down to 90.8%.

Examples

Y = exprnd(1, 40000, 5); r = sim_firquest(Y, 0.9, 0.05); [r.lower r.estimate r.upper] % brackets -log(0.1) = 2.3026

Reference: A. Lolos, C. Alexopoulos, D. Goldsman, K. D. Dingec, A. C. Mokashi, J. R. Wilson, “A Fixed-Sample-Size Procedure for Estimating Steady-State Quantiles Based on Independent Replications”, Proc. Winter Simulation Conference, 2025.

See also OA_FQUEST, OA_STS_QUANTILE_AREAS, OA_VONNEUMANN, OA_SHAPIROWILK