api.infer

infer_nhpp_ks(times, T, varargin)

INFER_NHPP_KS KS test for a non-homogeneous Poisson arrival process.

RESULT = INFER_NHPP_KS(TIMES, T) tests whether the arrival times TIMES in [0,T] came from an NHPP.

THE CONDITIONAL-UNIFORM TRANSFORMATION. Conditional on the number of arrivals in the interval, the arrival times of an NHPP are distributed as the order statistics of iid variables with cdf Lambda(t)/Lambda(T). Mapping the data through that cdf therefore turns ANY NHPP, whatever its rate, into iid uniforms, and one KS test then covers every rate function. Without a cumulative rate the rate is taken constant on the interval, which is the piecewise-constant approximation the reference uses on each subinterval.

WHY THE PLAIN TEST IS WEAK, AND WHAT FIXES IT. The CU KS test has “remarkably little power” against processes with non-exponential interarrival times, because it looks at the POSITIONS of the points and those stay nearly uniform for many non-Poisson processes. Lewis (1965) applies the Durbin (1961) transformation first: reorder the GAPS between the uniforms ascending, rescale each by how many gaps remain, and cumulate. That turns a difference in the gap DISTRIBUTION – exactly what a non-exponential renewal process has – into a difference in position, which KS can see. Measured on 400 replications of an Erlang-4 renewal process, the CU test rejects at its own size while the Lewis test rejects essentially always.

Options: ‘cumRate’ (the cumulative rate Lambda(t)), ‘method’ (‘lewis’ by default, ‘cu’ for the plain test), ‘T0’ (the left end, default 0).

Returns a struct with fields statistic, pvalue, n, uniforms and transformed.

Reference: S.-H. Kim, W. Whitt (2014). Are call center and hospital arrivals well modeled by nonhomogeneous Poisson processes? Manufacturing and Service Operations Management 16(3), 464-480; J. Durbin (1961), Biometrika 48, 41-55; P. A. W. Lewis (1965), JRSS B 27, 417-432.

See also TRACE_IDI, MAP_FIT.

mexify_infer

@brief MATLAB Coder script to generate MEX functions for infer_ module.

This script generates MEX (MATLAB Executable) versions of demand- inference functions for improved performance.

Skipped functions (Coder-incompatible):

infer_get_qlen_arrival - Cell-array data argument with ragged: per-class samples.
infer_gibbs - Cell-array data argument and use of: exist() for default-arg detection.
infer_rps - Calls lsqnonneg() (Optimization: Toolbox), not supported by Coder.

infer_mlps, infer_minps, infer_minps_setup, infer_fmlps - Take a Network model OOP object and

call its method getNumberOfServers().

infer_quick_model - Construct Network/Source/Queue/JobClass: OOP objects.
infer_fluid_ps_rt_likelihood - Returns a function handle (ode_h) and: consumes the sn struct.
sn_set_service_coc - Uses ProcessType enum, cell-of-cells: fields of the sn struct.

See also CODER, CODER.CONFIG, CODER.TYPEOF, CODEGEN.

infer_variational(spec, options)

OUT = INFER_VARIATIONAL(SPEC, OPTIONS)

Variational inference for Markovian queueing networks, following I. Perez, G. Casale, “Variational Inference for Markovian Queueing Networks”, Advances in Applied Probability 53(3), 2021.

The network trajectory is reparameterised by the transition counts Y^eta, eta=(i,j,c), so that the station marginals decouple:

x_{i,c}(t) = x_{i,c}(0) + sum_{eta in In(i,c)} Y^eta(t)

sum_{eta in Out(i,c)} Y^eta(t)

The variational family is a product of inhomogeneous pure-birth processes, one per transition, with rate nu^eta(t,y), times a product of Gamma densities over the unknown service rates. The state space is expanded by adding DELTA to every feasible rate, so that queue lengths may go negative and the approximating measure stays mutually absolutely continuous with the target (Sec. 3 of the paper); the original model is recovered as DELTA -> 0.

Each iteration performs, per transition, a backward pass for the Lagrange multipliers r^eta (Eq. 14) with multiplicative jumps at the observation epochs, the rate update

nu^eta(t,y) = exp(E log Xi^eta(t,y)) r^eta(t,y+1)/r^eta(t,y),

and a forward pass of the master equation for the marginal. The conjugate Gamma posteriors are then refreshed from the expected number of firings and the expected exposure time of each station-class pair.

Expectations over the other transitions are taken on a deterministic Halton lattice mapped through the inverse marginal c.d.f., so the estimator carries no random-number stream and is reproducible across the MATLAB, Java, Python and C++ implementations.

SPEC fields (plain arrays; (m,r) pairs are flattened column-major, so that pair (m,r) sits at index (r-1)*M+m):

arcs (T x 3) [i j c] transitions; i=0 external source, j=0 sink x0 (M x R) initial queue lengths sched (M x 1) 0=INF, 1=shared server (PS/FCFS), 2=external nservers (M x 1) number of servers routeprob (T x 1) routing probability p^c_{i,j} of each transition arcparam (T x 1) index in 1..P of the rate governing the arc, 0 if known arcrate (T x 1) known rate for arcs with arcparam==0, NaN otherwise alpha0 (P x 1) Gamma prior shapes beta0 (P x 1) Gamma prior rates obsTimes (K x 1) observation epochs obsData (K x M*R) observed queue lengths, NaN where not observed obsRange (M x R) support size of the uniform contamination epsilon scalar probability that a reading is faulty capacity (M x R) optional upper bound on the queue length, Inf by

default. In a closed network this is the chain population, and clamping the load there keeps the expanded state space from crediting a station with more jobs than the network holds

OPTIONS fields: tmax, dt, ngrid, ymax, nsamples, iter_max, tol, delta, floor, rate_max, verbose.

OUT fields: alpha, beta, rates (posterior mean rates), meanServiceTime, bound (ELBO per iteration), Y (T x G x ymax+1 marginals), nu, tgrid, qlen (G x M*R expected queue lengths), iter, converged, tailmass.

infer_mlps(model, node, rt, class, ql)

INFER_MLPS MLPS demand estimation using sn struct-level operations.

Estimates service demands at a PS queue using the Maximum Likelihood for Processor Sharing method. Pre-builds augmented CTMC models for each unique (tagClass, aQueue) combination once, then uses sn_set_service + solver_ctmc directly inside the optimization loop.

Inputs:: model - LINE Network model with delay rates set and queue rates to estimate node - PS queue node (Station object) rt - response time samples (column vector, n x 1) class - class of each sample (column vector, n x 1) ql - queue lengths at arrival (n x R matrix, per-class)

Returns:: demandEst - 1 x R vector of estimated mean service demands

infer_gibbs(data, nbCores, tol)

infer_lqn_setparams(model, paramSpec, a)

INFER_LQN_SETPARAMS Apply a parameter vector to a LayeredNetwork.

MODEL = INFER_LQN_SETPARAMS(MODEL, PARAMSPEC, A) sets the LQN parameters named in PARAMSPEC to the values in the vector A, in place, and returns the (same handle) MODEL. This is the parameter-injection half of the LQN parameter identification method (see INFER_LQN).

PARAMSPEC is a struct array; PARAMSPEC(i) has fields:

.type‘hostdem’ (activity mean host demand) or
‘think’ (task mean think time)

.name : name of the activity (‘hostdem’) or task (‘think’)

Mean values are injected as exponential distributions (SCV = 1), which is the assumption of the CASCON 2005 tracking method. The cached layered struct (MODEL.lsn) is invalidated so the next getStruct/solve re-reads the mutated model objects.

infer_lqn_optget(s, f, d)

INFER_LQN_OPTGET Return option field s.(f), or default d when absent/empty.

V = INFER_LQN_OPTGET(S, F, D) returns S.(F) if S is a struct with a non-empty field F, otherwise the default value D. Small helper shared by the LQN parameter identification routines (see INFER_LQN).

infer_lqn_jacobian(hfun, a, fdStep, fdFloor)

INFER_LQN_JACOBIAN Forward finite-difference sensitivity matrix of h at a.

[H, H0] = INFER_LQN_JACOBIAN(HFUN, A, FDSTEP, FDFLOOR) returns the sensitivity matrix H = dh/da evaluated at the parameter vector A by forward finite differences, and H0 = HFUN(A). HFUN maps a parameter vector to a numeric observation vector z = h(a). This is the approximate sensitivity matrix H_k used in the EKF update (Zheng, Yang, Woodside, Litoiu, Iszlai, “Tracking Time-Varying Parameters in Software Systems with Extended Kalman Filters”, CASCON 2005).

Each column is H(:,i) = (HFUN(a + d_i) - H0)/d_i with the per-parameter step d_i = FDSTEP * max(|a_i|, FDFLOOR). HFUN is evaluated numel(A)+1 times.

FDSTEP : relative perturbation (default 1e-3) FDFLOOR : minimum absolute perturbation scale (default 1e-6)

infer_lqn_getobs(names, metrics, obsSpec)

INFER_LQN_GETOBS Extract an observation vector from solved LQN metrics.

Z = INFER_LQN_GETOBS(NAMES, METRICS, OBSSPEC) returns the column vector Z of performance measures selected by OBSSPEC from the per-element numeric metric vectors in METRICS. This is the measurement-model half of the LQN parameter identification method (see INFER_LQN): it maps a solved model to the observation vector z = h(a).

NAMES : cell array of element names, as in LayeredNetworkStruct.names. METRICS : struct with fields QLen, Util, RespT, Tput (each a vector indexed

by element, aligned with NAMES, as returned by getEnsembleAvg).

OBSSPECstruct array; OBSSPEC(i) has fields:
.metric : ‘RespT’ | ‘Util’ | ‘Tput’ | ‘QLen’ .name : element name (processor for utilization; entry or reference

task for response time; entry/task for throughput)

infer_lqn_findbyname(container, name)

INFER_LQN_FINDBYNAME Locate an LQN element by name in a cell container.

EL = INFER_LQN_FINDBYNAME(CONTAINER, NAME) returns the first element of the cell array CONTAINER whose .name equals NAME, or [] if none matches. Used by the LQN parameter identification helpers to resolve activities and tasks referenced by name in a parameter specification (see INFER_LQN).

infer_lqn_ekf(hfun, a0, P0, Z, Q, R, options)

INFER_LQN_EKF Extended Kalman Filter for LQN parameter identification.

[AHAT, INFO] = INFER_LQN_EKF(HFUN, A0, P0, Z, Q, R, OPTIONS) tracks a hidden parameter vector across the measurement sequence Z using an Extended Kalman Filter, following Zheng, Yang, Woodside, Litoiu, Iszlai, “Tracking Time-Varying Parameters in Software Systems with Extended Kalman Filters”, CASCON 2005 (equations 1-9). The parameter is modelled as a zero-mean random walk a_k = a_{k-1} + w and the measurement as z_k = h(a_k) + v, where h is the (nonlinear) LQN performance model supplied as HFUN.

Inputs:

HFUNhandle mapping a parameter vector to a predicted observation
vector z = h(a) (evaluates the LQN model).

A0 : (np x 1) initial parameter estimate. P0 : (np x np) initial estimation-error covariance. Z : (no x nsteps) measurement matrix, one column per step. Q : (np x np) parameter-drift (process-noise) covariance. R : (no x no) measurement-error covariance. OPTIONS : struct with optional fields:

fdStep (1e-3), fdFloor (1e-6) - finite-difference steps clampPositive (true) - clamp estimates to > fdFloor aTrue ([]) - ground truth for Ea metric verbose (false)

Outputs:
AHAT : (np x nsteps) parameter estimate trajectory. INFO : struct with fields P (final covariance), Phist (per-step

covariances), e (no x nsteps prediction errors), zpred (predicted measurements), Er (prediction RMS), Ea (parameter tracking RMS vs aTrue when supplied, else []).

infer_lqn(model, paramSpec, obsSpec, Z, options)

INFER_LQN Identify hidden LQN parameters from measured performance data.

[MODEL, INFO] = INFER_LQN(MODEL, PARAMSPEC, OBSSPEC, Z, OPTIONS) estimates the LQN parameters named in PARAMSPEC (activity host demands and/or task think times) of the LayeredNetwork MODEL from the sequence of performance measurements Z, using an Extended Kalman Filter over the observation model defined by OBSSPEC. It implements Zheng, Yang, Woodside, Litoiu, Iszlai, “Tracking Time-Varying Parameters in Software Systems with Extended Kalman Filters”, CASCON 2005. A single measurement column with OPTIONS.QFac = 0 reduces to one-shot least-squares calibration.

Inputs:
PARAMSPEC : struct array of parameters to estimate; see INFER_LQN_SETPARAMS. OBSSPEC : struct array of observed metrics; see INFER_LQN_GETOBS. Z : (no x nsteps) measurements, no == numel(OBSSPEC). OPTIONS : struct with optional fields (defaults in parentheses):

solver (@SolverLN) observation-model solver constructor solveropts ([]) options struct passed to the solver QFac (0.1) drift-noise factor, Q_ii=(QFac*a0_i*cvA)^2 RFac (0.2) meas.-noise factor, R_ii=((RFac*zbar_i)/1.96)^2/gammaT cvA (1) parameter drift coefficient of variation gammaT ([]) T/Tstar ratio for R (else from T and Tstar, else 1) T, Tstar ([]) measurement interval and system constant Q, R ([]) explicit covariances (override the above) P0 ([]) initial covariance (default diag((0.5*a0)^2)) a0 ([]) initial estimate (default: current model values) aTrue ([]) ground truth, enables the Ea RMS metric fdStep (1e-3), fdFloor (1e-6), clampPositive (true), verbose (false)

Outputs:
MODEL : the LayeredNetwork with the final parameter estimate applied. INFO : EKF result struct (see INFER_LQN_EKF) augmented with ahat, a0,

Q, R, P0.

sn_set_service_coc(sn, stationIdx, classIdx, rate, scv)

SN_SET_SERVICE_COC Update service rate preserving cell-of-cells format.

Works like sn_set_service with autoRefresh=true, but writes mu, phi, proc, pie using cell-of-cells indexing (sn.mu{i}{k}) instead of 2D cell indexing (sn.mu{i,k}). This avoids reshaping the cell arrays, which would break solvers that expect cell-of-cells format.

Inputs:: sn - NetworkStruct stationIdx - station index (1-based) classIdx - class index (1-based) rate - new service rate (positive scalar) scv - squared coefficient of variation (default 1.0)

infer_rps(rt, class, ql, V)

INFER_RPS Regression for Processor Sharing (RPS) demand estimation

RT: response time samples (column vector with all samples) CLASS: class of each request sample (column vector) QL: queue length samples (matrix with R columns containing the

number of jobs of each class observed by each sample)

V: number of cores/servers

Based on mean-value analysis for PS stations. For a single server:

E[R_r] = E[D_r] * E[Q_bar_A]

where Q_bar_A is the total number of jobs seen upon admission, including the arriving job itself. For V > 1 processors, the queue length is split equally among the servers:

E[R_r] = E[D_r] * E[Q_bar_A] / V

The demand E[D_r] is estimated via non-negative least squares regression of response times against Q_bar_A / V.

Note: RPS assumes equal splitting among all V servers, which can lead to overestimation under low loads when fewer than V servers are actually busy. ERPS improves on this by using the mean number of busy servers instead of V.

infer_quick_model(open, stations, classes, servers, jobs, routing): QUICKMODEL Generate simple closed queueing network based on given parameters

infer_qmle(Q, N, Z)

infer_minps_setup(data, initSample, sampleSize, V, model, node)

MAIN_MINPS setups the input data for the MINPS estimation method and calls it

infer_minps(model, node, rt, class, ql)

INFER_MINPS MINPS demand estimation method.

Runs both MLPS and RPS estimators and selects the one with the smaller mean demand estimate.

Inputs:: model - LINE Network model with delay rates set and queue rates to estimate node - PS queue node (Station object) rt - response time samples (column vector, n x 1) class - class of each sample (column vector, n x 1) ql - queue lengths at arrival (n x R matrix, per-class)

Returns:: demandEst - 1 x R vector of estimated mean service demands

infer_get_qlen_arrival(data)

INFER_GET_QLEN_ARRIVAL Compute queue lengths at arrival from cell data.

Wrapper around infer_compute_ql_at_arrival for the legacy cell-based data format. Assumes data is available in standard format where data{3,k} contains arrival times (in ms) and data{4,k} contains response times for class k.

Inputs:: data - cell array in standard format (6 x K+1)

Returns:: ql - 1 x K cell array, each cell is numSamples(k) x K matrix

infer_fmlps(model, node, rt, class, ql, W)

INFER_FMLPS FMLPS demand estimation using sn struct-level operations.

Estimates service demands at a PS queue using the Fluid Maximum Likelihood for Processor Sharing method. Uses sn_set_service for fast parameter updates, avoiding model.reset()/getStruct() in the optimization loop.

Inputs:: model - LINE Network model with delay rates set and queue rates to estimate node - PS queue node (Station object) rt - response time samples (column vector, n x 1) class - class of each sample (column vector, n x 1) ql - queue lengths at arrival (n x R matrix, per-class) W - total population (number of threads/jobs)

Returns:: demandEst - 1 x R vector of estimated mean service demands

infer_fluid_ps_rt_likelihood(sn, taggedClass, y0_levels, Rsampled)

INFER_FLUID_PS_RT_LIKELIHOOD Fluid-based response time likelihood.

Builds an augmented model with tagged class K+1 by expanding the NetworkStruct arrays (following solver_fluid_passage_time’s pattern), then uses solver_fluid_odes to obtain the ODE handle.

Usage modes:

[ode_h, q_indices, augPhases] = infer_fluid_ps_rt_likelihood(sn, taggedClass): Build augmented model and return ODE handle + indexing info.
[~, ~, ~, LIKE] = infer_fluid_ps_rt_likelihood(sn, taggedClass, y0_levels, Rsampled): Build augmented model, solve ODE, and return likelihood.

Inputs:

sn - NetworkStruct (from model.getStruct()) taggedClass - class index of the tagged job y0_levels - M x K matrix of fluid levels per station per class (optional) Rsampled - observed response time (optional)

The augmented model adds class K+1 (tagged) with the same service rates as taggedClass. At the reference (queue) station, the tagged class absorbs by switching to taggedClass and routing to the delay station. The likelihood is extracted from the ODE derivative at the terminal state.

Reference: Casale, G. et al., “Fluid Analysis of Queueing in Processor Sharing Systems”

infer_compute_ql_at_arrival(at, at_jobid, rt, rt_jobid, class, R)

INFER_COMPUTE_QL_AT_ARRIVAL Compute per-class queue lengths at arrival.

Reconstructs the queue state seen by each arriving job using arrival and departure times. Arrival times and response times are matched by job ID, so they need not be in the same order or come from the same data source.

At ties, departures are processed before arrivals.

Inputs:

at - arrival times (column vector, n x 1) at_jobid - job IDs for arrival times (column vector, n x 1) rt - response times (column vector, m x 1, m >= n) rt_jobid - job IDs for response times (column vector, m x 1) class - class of each arrival sample (column vector, n x 1,

ordered consistently with at and at_jobid)

R - number of classes

Returns:

ql - n x R matrix of per-class queue lengths at each arrival,: rows ordered consistently with the input at/at_jobid