lang.processes 

An astract class for Markovian distributions

Method Summary

D(i, wantSparse): Di = D(i)

evalCDF(t): FT = EVALCDF(SELF,T)

evalLST(s): L = EVALLAPLACETRANSFORM(S)

getInitProb(): ALPHA = GETINITPROB()

getMean(): MEAN = GETMEAN()

getMoments(n): EXN = GETMOMENTS(N)

getMu(): MU = GETMU()

getNumberOfPhases(): PHASES = GETNUMBEROFPHASES()

getPhi(): PHI = GETPHI()

getProcess(): D = GETPROCESS()

getRepresentation(): D = GETREPRESENTATION() Return the (D0,D1,…) matrix representation of the process. Alias of getProcess, mirroring the Java/Python API.

getSCV(): SCV = GETSCV() Get the squared coefficient of variation of the distribution (SCV = variance / mean^2)

getSkewness(): SKEW = GETSKEWNESS()

getSubgenerator(): T = GETSUBGENERATOR()

isImmediate(): BOOL = ISIMMEDIATE() Check if the distribution is equivalent to an Immediate distribution

plot(): PLOT()

sample(n): X = SAMPLE(N)

setMean(MEAN): SELF = UPDATEMEAN(MEAN)

setProcess(D)

setRate(RATE): SETRATE(SELF,RATE)

class Trace

Bases: Replayer

Empirical time series from a trace, alias for Replayer

Constructor Summary

Trace(data): SELF = TRACE(data)

Method Summary

getFeatureName()

FEATNAME = GETFEATURENAME() ‘Trace’, not the inherited ‘Replayer’.

NAME stays ‘Replayer’ because it selects the ProcessType and the JSON wire type, and a Trace replays samples exactly as a Replayer does. The registry name is what makes the Trace entry reachable; a solver declaring only ‘Replayer’ still accepts the model, through SolverFeatureSet.generalizationOf.

getMoments()

class MarkovProcess

Bases: Process

A class for a continuous time Markov chain

Constructor Summary

MarkovProcess(InfGen, isFinite, stateSpace): SELF = MARKOVPROCESS(InfGen, isInfinite, stateSpace)

Property Summary

infGen

isfinite

stateSpace

Method Summary

aggregate(MS, method, param): [P, EPS, EPSMAX] = AGGREGATE(MS, METHOD, PARAM) Aggregation-disaggregation over the macrostate partition MS, a cell array of state-index vectors. METHOD is ‘courtois’ (PARAM is the randomization rate q), ‘kms’ or ‘takahashi’ (PARAM is the iteration count, default 10), or ‘multi’ (PARAM is the second-level partition MSS). EPS is the nearly-complete- decomposability index of the partition and EPSMAX the largest index for which the approximation is meant to hold.

static fromSampleSysAggr(sa)

getGenerator(): Q = GETGENERATOR()

getProbState(state): Use Cramer’s rule to compute the probability of a single state

hittingTime(targetStates): H = HITTINGTIME(TARGETSTATES) Mean TIME to reach any state in TARGETSTATES, zero on the target set itself and Inf from a state that cannot reach it. The twin MarkovChain.hittingTime counts STEPS instead, so the two answer different questions about the same jump structure.

isFeasible(): BOOL = ISFEASIBLE() True when the generator is a valid one.

plot()

plot3()

static rand(nStates): creates a random CTMC

sample(n)

sens(dQ): DPI = SENS(DQ) Sensitivity of the stationary distribution to a scalar parameter, given the derivative DQ of the generator.

setStateSpace(stateSpace)

solve()

solveRelative(refstate): P = SOLVERELATIVE(REFSTATE) Equilibrium distribution relative to REFSTATE, i.e. with p(REFSTATE)=1. Unnormalized by construction, so it is defined even where the normalizing constant is not.

stochComp(I): S = STOCHCOMP(I) Stochastic complement of the states I, a generator on that subset. Use stochCompFull to also obtain the partitioned blocks.

stochCompFull(I): [S, Q11, Q12, Q21, Q22, T] = STOCHCOMPFULL(I) Stochastic complement of the states I together with the blocks of the generator partitioned by I and its complement, and the return-path term T = Q12*inv(-Q22)*Q21, so that S = Q11 + T.

timeAverage(pi0, t): [PITIMEAVG, PIEXIT] = TIMEAVERAGE(PI0, T) Time-averaged distribution over [0,T] and its endpoint.

toDTMC(q): TODTMC - Alias for toMarkovChain for backwards compatibility

toEmbedded(): A = TOEMBEDDED() Embedded jump chain, i.e. the DTMC of the states visited at transition epochs. Unlike toDTMC (uniformization) it does not preserve the stationary distribution, since it drops the holding times; an absorbing state stays absorbing.

toMarkovChain(q)

toTimeReversed()

transient(pi0, t, method): [PI_T, KMAX] = TRANSIENT(PI0, T, METHOD) Distribution at time T from PI0 (uniform if empty). METHOD is ‘unif’ (Jensen uniformization, the default) or ‘foxglynn’, whose weights avoid evaluating the Poisson terms directly.

transientProb(pi0, t): [PI_T, KMAX] = TRANSIENTPROB(PI0, T) Alias of transient, under the name the JAR must use since ‘transient’ is a Java keyword.

class Expolynomial

Expolynomial distribution with density f(x) = sum ci * x^ai * exp(-li*x)

Represents an expolynomial density over a bounded domain [eft, lft], matching the GEN expolynomial format of external stochastic Petri net tools.

Constructor Summary

Expolynomial(density, eft, lft)

EXPOLYNOMIAL Create an Expolynomial distribution instance

@param density Density expression string in expolynomial (GEN) format @param eft Earliest firing time (lower bound of support) @param lft Latest firing time (upper bound of support, can be Inf) @return self Expolynomial distribution instance

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T)

getMean(): EX = GETMEAN()

getProcess(): PROC = GETPROCESS()

getSCV(): SCV = GETSCV()

sample(n): X = SAMPLE(N)

class Distribution

Bases: Copyable

Distribution Abstract base class for statistical distributions

Distribution provides the common interface and functionality for all statistical distributions used in queueing models. It defines abstract methods for sampling, computing moments, and obtaining distribution properties that must be implemented by concrete distribution classes.

@brief Abstract base class for all statistical distribution types

Key characteristics: - Abstract interface for all distributions - Support for both discrete and continuous distributions - Moment computation (mean, variance, SCV) - Random sampling capabilities - Parameter management and validation - Copyable interface for cloning distributions

Distribution hierarchy includes: - Continuous distributions (Exp, Erlang, HyperExp, etc.) - Discrete distributions (Poisson, Bernoulli, Zipf, etc.) - Phase-type distributions (PH, APH, Coxian, etc.) - Markov-modulated processes (MAP, MMPP, etc.)

Common usage patterns: @code dist = Exp(1.5); % Exponential with rate 1.5 mean_val = dist.getMean(); % Get mean value samples = dist.sample(100); % Generate 100 samples @endcode

Property Summary

name

obj

params

support: support interval

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T) Evaluate the cumulative distribution function at t

evalLST(s): L = EVALLAPLACETRANSFORM(S) Evaluate the Laplace transform of the distribution function at t

evalProbInterval(t0, t1): DELTA = EVALPROBINTERVAL(SELF,T0,T1) Evaluate the probability mass between t0 and t1 (t1>t0)

getFeatureName()

FEATNAME = GETFEATURENAME() The SolverFeatureSet entry this distribution is marked under by getUsedLangFeatures.

Separate from NAME because that one also selects the ProcessType and the JSON wire type: a subclass whose registry name is more specific than its process type (a Trace, or a two-phase Coxian) says so here without moving to a different process type. Defaults to NAME, so a distribution needing no distinction is unaffected.

getMean(): MEAN = GETMEAN() Get distribution mean

getNumParams(): NPARAM = GETNUMPARAMS() Returns the number of parameters needed to specify the distribution

getParam(id): PARAM = GETPARAM(SELF,ID) Return the parameter associated to the given id

getRate()

GETRATE Get the rate parameter (inverse of mean)

@brief Returns the rate parameter defined as 1/mean @return RATE Rate parameter (1/mean)

getSCV(): SCV = GETSCV() Get distribution squared coefficient of variation (SCV = variance / mean^2)

getSkewness(): SKEW = GETSKEWNESS() Get distribution skewness

getVar(): VAR = GETVAR() Get distribution variance

isContinuous(): BOOL = ISCONTINUOUS() Check if the distribution is discrete

isDisabled()

BOOL = ISDISABLED(): Check if the distribution is equivalent to a Disabled distribution

bool = cellfun(@(c) isnan(c.paramValue), self.params)

isDiscrete(): BOOL = ISDISCRETE() Check if the distribution is discrete

isImmediate(): BOOL = ISIMMEDIATE() Check if the distribution is equivalent to an Immediate distribution

sample(n)

SAMPLE Generate random samples from the distribution

@brief Generates n random samples from this distribution @param n Number of samples to generate @return X Vector of n random samples

setNumParams(numParam): NPARAM = SETNUMPARAMS(NUMPARAM) Initializes the parameters

setParam(id, name, value): SETPARAM(ID, NAME, VALUE, TYPECLASS) Set a distribution parameter given id, name, value, Java class type (for JMT translation)

class DiscreteSampler

A class for discrete distributions specified from the probability mass function

Constructor Summary

DiscreteSampler(p, x): SELF = DISCRETESAMPLER(P,X)

Method Summary

evalCDF(k): FT = EVALCDF(SELF,K)

evalPMF(v): PK = EVALPMF(V)

getMean(): EX = GETMEAN()

getSCV(): SCV = GETSCV()

isDisabled(): BOOL = ISDISABLED()

sample(n): X = SAMPLE(N)

class Coxian

Bases: Markovian

The coxian statistical distribution

Constructor Summary

Coxian(varargin): SELF = COXIAN(VARARGIN) Constructs a Coxian distribution from phase rates and completion probabilities, with entry probability 1 on the first phase

Method Summary

static fitCentral(MEAN, VAR, SKEW): CX = FITCENTRAL(MEAN, VAR, SKEW)

static fitMeanAndSCV(MEAN, SCV): [CX,MU,PHI] = FITMEANANDSCV(MEAN, SCV) Fit a Coxian distribution with given mean and squared coefficient of variation (SCV=variance/mean^2)

getFeatureName()

FEATNAME = GETFEATURENAME() ‘Cox2’ at two phases, ‘Coxian’ otherwise.

The registry carries both names, and the Cox2 entry can only mean the two-phase Coxian here: Cox2 is a static factory returning a Coxian, so there is no Cox2 object to mark, and without this the entry is unreachable. The C++ port already types every two-phase Coxian as ProcessType.COX2, so the reading is shared. A solver declaring just ‘Coxian’ stays accepting through SolverFeatureSet.generalizationOf.

getMean(): EX = GETMEAN() Get distribution mean

getMu(): MU = GETMU() Get vector of rates

getNumberOfPhases(): PHASES = GETNUMBEROFPHASES() Return number of phases in the distribution

getPhi(): PHI = GETPHI() Get vector of completion probabilities

getSCV(): SCV = GETSCV() Get the squared coefficient of variation of the distribution (SCV = variance / mean^2)

class Uniform

Uniform Continuous uniform distribution over a bounded interval

Uniform represents the continuous uniform distribution with equal probability density over a specified interval [min, max]. This distribution models scenarios where all values within a range are equally likely, making it useful for modeling bounded random processes with no preferred values.

@brief Uniform distribution with equal probability over [min, max] interval

Key characteristics: - Two parameters: minimum (min) and maximum (max) values - Mean = (min + max) / 2 - Variance = (max - min)² / 12 - Constant probability density over [min, max] - Support: [min, max]

The uniform distribution is used for: - Modeling bounded processes with no bias - Random number generation in simulations - Representing equally likely outcomes - Initial parameter estimation - Load testing with uniform traffic patterns

Example: @code service_dist = Uniform(1.0, 3.0); % Service time between 1 and 3 units % Mean = 2.0, all values in [1.0, 3.0] equally likely samples = service_dist.sample(1000); @endcode

Constructor Summary

Uniform(minVal, maxVal)

UNIFORM Create a Uniform distribution instance

@brief Creates a Uniform distribution over the interval [minVal, maxVal] @param minVal Minimum value of the distribution support @param maxVal Maximum value of the distribution support (must be > minVal) @return self Uniform distribution instance

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T)

evalLST(s): L = EVALST(S)

getMean(): EX = GETMEAN()

getProcess(): PROC = GETPROCESS()

getSCV(): SCV = GETSCV()

sample(n): X = SAMPLE(N)

class PHt

PHt Time-inhomogeneous phase-type distribution (Ph_t).

Following Ko and Pender (Oper. Res. Lett. 45, 2017), a Ph_t is an ordinary phase-type distribution whose initial vector and sub-generator are functions of the wall clock, alpha(t) and S(t), required only to be locally integrable. This class realises that definition with a piecewise-constant schedule: segment k covers [breakpoints(k), breakpoints(k+1)) and carries the pair (alpha{k}, S{k}), so breakpoints has one more entry than the lists. The exit vector is s(t) = -S(t)e.

Because both the phase and the elapsed service depend on absolute time, a Ph_t service time is a function of the epoch at which service starts: sampleFrom(t0) is the operative sampler, and sample walks one path.

Setting h = 1 with S = -mu_k recovers a time-varying exponential, whose completion stream at a saturated server is the NHPP with rates mu_k.

Like MAPt this does NOT extend Markovian, so that isMarkovian-gated code cannot read it as a single stationary (alpha, S) pair; and the scalar summaries getSCV, getSkewness, evalCDF, evalLST return NaN, the distribution of a service time being different at every start epoch.

Constant support. The fluid solver expresses a segment as a per-entry multiplier on the time-averaged nominal, so the constructor requires the sparsity pattern of alpha, of the off-diagonal S and of the exit vector to be identical across segments.

The process representation stores:: process{1} : 1-by-(n+1) row vector of breakpoints process{2} : 1-by-n cell of alpha row vectors process{3} : 1-by-n cell of S matrices process{4} : logical, true if cyclic

Constructor Summary

PHt(breakpoints, alpha, S, cyclic): SELF = PHT(BREAKPOINTS, ALPHA, S, CYCLIC)

Property Summary

S: 1-by-n cell of sub-generators

alpha: 1-by-n cell of initial probability row vectors

breakpoints: 1-by-(n+1) segment boundaries, strictly increasing

cyclic: logical, whether the schedule repeats

process: {breakpoints, alpha, S, cyclic}, for serialization

sampleClock: wall-clock position of the next sample (see sample)

Method Summary

evalCDF(t): F = EVALCDF(T) NaN; see getSCV.

evalLST(s): L = EVALLST(S) NaN; see getSCV.

getAlphaAt(t): A = GETALPHAAT(T) alpha in force at T; the last segment’s vector past a non-cyclic horizon, where getSAt is zero anyway.

getAlphaSegments()

getBreakpoints()

getMean(): MEAN = GETMEAN() Mean of the time-averaged phase-type.

getNumSegments()

getNumberOfPhases()

getPeriod(): T = GETPERIOD() Horizon length, which is the period when cyclic.

getProcess()

getRateAt(t): R = GETRATEAT(T) Completion rate of the PH in force at T.

getRateSchedule(): SCHED = GETRATESCHEDULE() The parameterisation of the process; the scalar interval summaries are not. Model compilation recognises a schedule-bearing process by this method rather than by class name.

getSAt(t): M = GETSAT(T) S in force at T; zero past a non-cyclic horizon.

getSCV(): SCV = GETSCV() NaN: the service-time distribution differs at every start epoch, so there is no single i.i.d. law for an SCV to summarise. Returning the SCV of the time-averaged representation would report a time-varying process as a stationary one to every consumer of sn.scv.

getSSegments()

getSegmentIndexAt(t): IDX = GETSEGMENTINDEXAT(T) Active segment, 0 past a non-cyclic horizon.

getSkewness(): SKEW = GETSKEWNESS() NaN; see getSCV.

getTimeAverageProcess(): [ABAR, SBAR] = GETTIMEAVERAGEPROCESS() Width-weighted average over the horizon. A convex combination of sub-generators is a sub-generator and of probability vectors a probability vector, so the nominal is a valid phase-type.

getTimeAverageProcessMAP(): [D0BAR, D1BAR] = GETTIMEAVERAGEPROCESSMAP() The nominal as a (D0, D1) pair, D1 = s*alpha, for the fluid carrier.

getTimeAverageRate(): RATE = GETTIMEAVERAGERATE() Completion rate of the time-averaged PH.

isCyclic()

resetSampleClock(): RESETSAMPLECLOCK() Restart the sample path at the schedule start.

sample(n): X = SAMPLE(N) Draws N successive service times along ONE sample path: sample i starts where sample i-1 completed, not at a fixed epoch. Use resetSampleClock to restart, or sampleFrom to draw a service time starting at a chosen epoch.

sampleFrom(t0): X = SAMPLEFROM(T0) Service time for a job whose service starts at wall clock T0. Exact: within a segment the phase process is a homogeneous absorbing CTMC, and by the memoryless property the residual holding time may be redrawn at a breakpoint. The initial phase is drawn from alpha in force at T0. Returns 0 when a non-cyclic horizon is exhausted before absorption.

class MarkovChain

Bases: Process

An abstract class for a discrete time Markov chain

Constructor Summary

MarkovChain(transMat, isFinite): SELF = MARKOVCHAIN(transMat, isInfinite)

Property Summary

isfinite

stateSpace

transMat

Method Summary

static fromSampleSysAggr(sa)

getTransMat()

hittingTime(targetStates): H = HITTINGTIME(TARGETSTATES) Mean number of steps to reach any state in TARGETSTATES.

isFeasible(): BOOL = ISFEASIBLE() True when the transition matrix is stochastic.

plot()

static rand(nStates): creates a random DTMC

sample(n): STS = SAMPLE(N) - Simulate n steps of the DTMC from a random initial state

setStateSpace(stateSpace)

solve(): PI = SOLVE() Stationary distribution of the DTMC. Twin of MarkovProcess.solve.

stochComp(I): S = STOCHCOMP(I) Stochastic complement of the states I, a DTMC on that subset. Use stochCompFull to also obtain the partitioned blocks.

stochCompFull(I): [S, P11, P12, P21, P22] = STOCHCOMPFULL(I) Stochastic complement of the states I together with the blocks of the transition matrix partitioned by I and its complement.

toCTMC(): TOCTMC - Alias for toMarkovProcess for backwards compatibility

toMarkovProcess()

toTimeReversed()

transient(pi0, steps): PI_T = TRANSIENT(PI0, STEPS) Distribution at each step 0,…,STEPS from PI0 (uniform if empty).

transientProb(pi0, steps): PI_T = TRANSIENTPROB(PI0, STEPS) Alias of transient, under the name the JAR must use since ‘transient’ is a Java keyword.

transientUnif(pi0, t): [PI_T, KMAX] = TRANSIENTUNIF(PI0, T) Distribution at time T of the DTMC seen through uniformization. The chain is read as the randomized image of a CTMC, so T is continuous here, unlike the step count taken by transient.

class MAPt

MAPt Time-inhomogeneous Markovian arrival process (MAP_t).

Following Ko and Pender (Oper. Res. Lett. 45, 2017), a MAP_t is an ordinary MAP whose two matrices are functions of the wall clock, D0(t) and D1(t), required only to be locally integrable. This class realises that definition with a piecewise-constant schedule, which is dense in L1_loc and is the form that serialises: segment k covers [breakpoints(k), breakpoints(k+1)) and carries the pair (D0{k}, D1{k}), so breakpoints has one more entry than the matrix lists. D0 holds transition rates without an arrival, D1 the rates that generate one, and D0+D1 is a generator in every segment.

Two horizon conventions, as for NHPP:

cyclicthe schedule repeats with period: T = breakpoints(end) - breakpoints(1).
non-cyclicoutside [breakpoints(1), breakpoints(end)) the process is: frozen in its last phase and emits nothing, so a non-cyclic MAP_t is a transient construct.

Setting h = 1 with D0 = -lambda_k, D1 = lambda_k recovers exactly the NHPP with the same breakpoints and rates.

This is neither a renewal process nor a time-homogeneous one, so the scalar summaries that presuppose an i.i.d. interval distribution – getSCV, getSkewness, evalCDF, evalLST – are undefined and return NaN rather than a representative value that would misreport the process as stationary. The schedule is the parameterisation: read it with getRateSchedule.

The class deliberately does NOT extend Markovian. Code gated on isMarkovian reads getProcess as a single stationary (D0, D1) pair and would silently drop the schedule; NHPP avoids the base class for the same reason.

Constant support. The fluid solver expresses a segment as a per-entry multiplier on the time-averaged nominal, so the constructor requires the sparsity pattern of the matrices to be identical across segments. A schedule that switches a transition on or off is refused outright rather than silently losing it.

The process representation stores:: process{1} : 1-by-(n+1) row vector of breakpoints process{2} : 1-by-n cell of D0 matrices process{3} : 1-by-n cell of D1 matrices process{4} : logical, true if cyclic

Solver support. SolverFLD honours the schedule in getTranAvg, and its ‘kp’ method integrates the Ko-Pender fluid and diffusion limits. Every other solver rejects a model using it via the standard unsupported-feature check.

Constructor Summary

MAPt(breakpoints, D0, D1, cyclic): SELF = MAPT(BREAKPOINTS, D0, D1, CYCLIC)

Property Summary

D0: 1-by-n cell of no-arrival rate matrices

D1: 1-by-n cell of arrival-generating rate matrices

breakpoints: 1-by-(n+1) segment boundaries, strictly increasing

cyclic: logical, whether the schedule repeats

process: {breakpoints, D0, D1, cyclic}, for serialization

sampleClock: wall-clock position of the next sample (see sample)

samplePhase: phase of the modulating chain at sampleClock

Method Summary

static checkCommonSupport(mats, ignoreDiagonal, label): CHECKCOMMONSUPPORT(MATS, IGNOREDIAGONAL, LABEL) Reject a schedule whose matrices do not share one sparsity pattern. The fluid solver expresses a segment as a per-entry multiplier on a nominal matrix, and that multiplier is undefined where the nominal entry is zero.

evalCDF(t): F = EVALCDF(T) NaN; see getSCV.

evalLST(s): L = EVALLST(S) NaN; see getSCV.

getBreakpoints()

getD0At(t): M = GETD0AT(T) D0 in force at T; zero past a non-cyclic horizon.

getD0Segments()

getD1At(t): M = GETD1AT(T) D1 in force at T; zero past a non-cyclic horizon.

getD1Segments()

getMean(): MEAN = GETMEAN() Palm mean interval of the time-averaged MAP.

getNumSegments()

getNumberOfPhases()

getPeriod(): T = GETPERIOD() Horizon length, which is the period when cyclic.

getProcess()

getRateAt(t): R = GETRATEAT(T) Arrival rate of the MAP in force at T, zero past a non-cyclic horizon. This is the stationary rate of that segment, not the instantaneous conditional intensity, which depends on the phase.

getRateSchedule(): SCHED = GETRATESCHEDULE() The parameterisation of the process; the scalar interval summaries are not. Model compilation recognises a schedule-bearing process by this method rather than by class name.

getSCV(): SCV = GETSCV() NaN: a MAP_t is neither renewal nor time-homogeneous, so there is no i.i.d. interval distribution for an SCV to summarise. Returning the SCV of the time-averaged MAP would report a time-varying process as a stationary one to every consumer of sn.scv.

getSegmentIndexAt(t): IDX = GETSEGMENTINDEXAT(T) Active segment, 0 past a non-cyclic horizon.

getSkewness(): SKEW = GETSKEWNESS() NaN; see getSCV.

getTimeAverageProcess(): [D0BAR, D1BAR] = GETTIMEAVERAGEPROCESS() Width-weighted average pair over the horizon. This is the stationary carrier of the phase structure where a solver needs a time-homogeneous one; a convex combination of generators is a generator, so it is itself a valid MAP.

getTimeAverageRate(): RATE = GETTIMEAVERAGERATE() Arrival rate of the time-averaged MAP. For h = 1 this is exactly the NHPP width-weighted average intensity.

isCyclic()

nextArrival(from, phase): [X, PHASE] = NEXTARRIVAL(FROM, PHASE) Time to the next arrival from wall clock FROM in the given phase, and the phase after that arrival. Exact: within a segment the phase process is a homogeneous CTMC, and by the memoryless property the residual holding time may be redrawn at a breakpoint, so the boundary is crossed by advancing the clock and resampling under the new matrices. Returns 0 when a non-cyclic horizon is exhausted.

resetSampleClock(): RESETSAMPLECLOCK() Restart the sample path at the schedule start.

sample(n): X = SAMPLE(N) Draws N successive interarrival times along ONE sample path. Both the intensity and the phase depend on absolute time, so this advances an internal clock and phase across calls. Use resetSampleClock to restart. A non-cyclic schedule that runs out returns 0 for every remaining sample.

class DiscreteUniform

The uniform statistical distribution

Constructor Summary

DiscreteUniform(minVal, maxVal): SELF = UNIFORM(MINVAL, MAXVAL)

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T)

evalLST(s): L = EVALST(S) Evaluate the Laplace-Stieltjes transform of the distribution function at s For DiscreteUniform(a, b), LST(s) = (e^(-as) - e^(-(b+1)s)) / ((b-a+1)(1 - e^(-s)))

evalPMF(k): P = EVALPMF(SELF,K)

getMean(): EX = GETMEAN()

getProcess(): PROC = GETPROCESS()

getSCV(): SCV = GETSCV()

sample(n): X = SAMPLE(N)

class DMAP

Discrete-time Markovian Arrival Process

Models discrete-time arrival streams with D0 and D1 sub-stochastic matrices. D0+D1 is a stochastic matrix (row sums = 1).

Constructor Summary

DMAP(D0, D1)

Method Summary

D(i, wantSparse)

evalMeanT(t)

getMean()

getNumberOfPhases()

getRate()

getSCV(): SCV = GETSCV()

getVar()

VAR = GETVAR() Variance of the inter-arrival count: 2*al*(I-D0)^-2*e - m - m^2.

DECLARED HERE because the inherited Markovian.getVar/getSCV call map_scv({D0,D1}), a CONTINUOUS-time formula that reads D0+D1 as a generator. For a DMAP that matrix is STOCHASTIC, so the stationary solve behind it is singular and the number it returned was not the variance of anything. jline.lang.processes.DMAP and the native Python DMAP already carry this formula; MATLAB was the outlier.

isImmediate()

static rand(order)

sample(n)

setMean(MEAN)

dist_scale_rate(distrib, factor)

Rate-scaled copy of a distribution, preserving its shape.

Returns a new distribution whose rate is FACTOR times the rate of DISTRIB, i.e. the time-scaled variable X/FACTOR. The scaling is exact: every moment of order n is divided by FACTOR^n, so the mean is divided by FACTOR while the SCV, the skewness and the whole shape of the distribution are preserved.

The scaled object is rebuilt from the parameters of the original, rather than by rescaling its Markovian representation, so that the parameter list stays coherent with the distribution family. Solvers that serialize the model (JMT, LDES) read the parameters, and would otherwise export the unscaled process.

This is the perturbation primitive of the finite-difference branch of getSensitivityTable: scaling the rate at a station-class by (1+h) is exactly the perturbation the derivative d(.)/d(rate) is taken along.

Parameters:

distrib – A Distribution object.
factor – Positive scaling factor for the rate.

Returns:

scaled – A new Distribution of the same family with rate*FACTOR.

dist_fit_me(mean, scv, maxPhases)

Fits a matrix exponential to a mean and squared coefficient of variation

For SCV < 1 the fit is the convolution X = c*Y + Z of a scaled concentrated matrix exponential Y (unit mean, minimal SCV sY for its order) with an independent exponential Z. Writing c + d = mean and c^2*sY + d^2 = scv*mean^2,

c = mean*(1 - sqrt(1 - (1+sY)*(1-scv)))/(1 + sY), d = mean - c,

so every target in [sY/(1+sY), 1] is matched EXACTLY in 2n+2 phases. The exponential tail is what makes the convolution reach up to SCV 1; the concentrated part is what makes it reach far below the Erlang bound 1/order at the same order.

The order is the smallest tabulated one that reaches the target, capped by maxPhases when given: with a phase budget an Erlang can only reach 1/maxPhases, while this construction reaches O(1/maxPhases^2), and the residual SCV is then the closest achievable from below.

SCV >= 1 is outside the range of a concentrated ME (its SCV never exceeds 0.34), and the caller keeps its own hyperexponential fit there.

Mirrors native Python fit_me_mean_scv and jline.lang.processes.MEFit.

Parameters:

mean – Target mean, positive
scv – Target squared coefficient of variation, in (0,1)
maxPhases – Optional cap on the number of phases, 0 for no cap

Returns:

me – ME distribution with the requested mean and, budget permitting, SCV

Examples

me = dist_fit_me(mean, scv)
me = dist_fit_me(mean, scv, maxPhases)

class RAP

Bases: Markovian

Rational Arrival Process (RAP) distribution

RAP is a generalization of the Markovian Arrival Process (MAP) where the matrices H0 and H1 represent hidden and visible transitions respectively, but with relaxed constraints compared to MAP.

Representation: - H0: matrix of hidden transition rates (transitions without arrivals) - H1: matrix of visible transition rates (transitions with arrivals) - H0 + H1 must form a valid infinitesimal generator (row sums = 0) - All eigenvalues of H0 must have negative real parts - Dominant eigenvalue of H0 must be negative and real

The marginal distribution of inter-arrival times is a Matrix Exponential (ME).

Constructor Summary

RAP(H0, H1)

RAP Create a Rational Arrival Process instance

@brief Creates a RAP with the given H0 and H1 matrices @param H0 Hidden transition matrix (square matrix) @param H1 Visible transition matrix (square matrix, same size as H0) @return self RAP distribution instance

Property Summary

H0: Hidden transition matrix

H1: Visible transition matrix

Method Summary

evalCDF(t)

FT = EVALCDF(SELF, T) Evaluate the cumulative distribution function at t

For RAP, the marginal CDF is same as MAP

evalLST(s): L = EVALST(S) Evaluate the Laplace-Stieltjes transform at s

evalPDF(t)

FT = EVALPDF(SELF, T) Evaluate the probability density function at t

For RAP, the marginal PDF is same as MAP

static fitMoments(moms)

RAP = FITMOMENTS(MOMS) Create RAP by fitting the given moments Uses BuTools RAPFromMoments algorithm

@param moms Array of moments @return rap RAP distribution matching the given moments

static fitMomentsAndCorrelations(moms, corrs)

RAP = FITMOMENTSANDCORRELATIONS(MOMS, CORRS) Create RAP by fitting the given moments and correlations Uses BuTools RAPFromMomentsAndCorrelations algorithm

@param moms Array of moments @param corrs Array of lag-k correlations @return rap RAP distribution matching the given moments and correlations

static fromErlang(k, rate)

RAP = FROMERLANG(K, RATE) Create RAP from Erlang renewal process Convenience method showing that Erlang is a special case of RAP

@param k Number of phases @param rate Rate parameter for each phase @return rap RAP distribution equivalent to Erlang(k, rate)

static fromMAP(map)

RAP = FROMMAP(MAP) Create RAP from Markovian Arrival Process Convenience method showing that MAP is a special case of RAP

@param map MAP distribution instance @return rap RAP distribution equivalent to the given MAP

static fromPoisson(rate)

RAP = FROMPOISSON(RATE) Create RAP from exponential renewal process (Poisson) Convenience method showing that Poisson is a special case of RAP

@param rate Arrival rate (lambda) @return rap RAP distribution equivalent to Poisson(rate)

getACF(lags): ACF = GETACF(SELF, LAGS) Get autocorrelation function at specified lags

getH0(): H0_OUT = GETH0() Get H0 matrix (hidden transitions)

getH1(): H1_OUT = GETH1() Get H1 matrix (visible transitions)

getIDC(t): IDC = GETIDC(SELF, T) Get index of dispersion for counts If t is not provided, returns asymptotic IDC

getMean(): MEAN_VAL = GETMEAN() Get mean of the RAP distribution

getNumberOfPhases(): PHASES = GETNUMBEROFPHASES() Get number of phases in the RAP representation

getProcess(): PROC = GETPROCESS() Get process representation {H0, H1}

getRate(): LAM = GETRATE() Get arrival rate (lambda) of the RAP

getSCV(): SCV = GETSCV() Get squared coefficient of variation

getVar(): VAR_VAL = GETVAR() Get variance of the RAP distribution

sample(n): X = SAMPLE(N) Get n samples from the distribution using rap_sample

class Prior

Bases: Distribution

Prior Discrete prior distribution over alternative distributions

Prior represents parameter uncertainty by specifying a discrete set of alternative distributions with associated probabilities. When used with setService or setArrival, it causes the UQ solver to expand the model into a family of networks, one for each alternative.

This is NOT a mixture distribution - each alternative represents a separate model realization with its associated prior probability.

Two forms are supported: - Discrete: an explicit set of alternative distributions with weights. - Continuous: a density f(theta) over a scalar parameter theta plus a

factory mapping theta to a Distribution. This is the form required by the epistemic uncertainty propagation of Trivedi and Bobbio (2017), Sec. 3.4, where the unconditional measure is the integral of the conditional measure against f(theta). The continuous form is reduced to a weighted alternative set by discretize(), so both forms are consumed identically downstream.

@brief Discrete or continuous prior for parameter uncertainty modeling

Key characteristics: - Discrete set of alternative distributions, or a continuous parameter density - Probability-weighted alternatives (must sum to 1) - Used with UQ solver for Bayesian analysis

Example: @code % Discrete form: service time with uncertain rate prior = Prior({Exp(1.0), Exp(2.0), Erlang(2,1.5)}, [0.4, 0.35, 0.25]); queue.setService(class, prior);

% Continuous form: rate is itself Erlang-distributed prior = Prior(Erlang(10, 3), @(lambda) Exp(lambda));

% Continuous form from k lifetime observations summing to s % (Jeffreys posterior of Trivedi-Bobbio Eq. 3.71) prior = Prior.fromSample(10, 5.0);

% Solve with UQ wrapper post = UQ(model, @SolverMVA); avgTable = post.getAvgTable(); % Prior-weighted expectations postTable = post.getPosteriorTable(); % Per-alternative breakdown @endcode

Constructor Summary

Prior(varargin)

PRIOR Create a prior distribution instance

@brief Creates a Prior in either discrete or continuous form

PRIOR(DISTRIBUTIONS, PROBABILITIES) discrete form. @param distributions Cell array of Distribution objects @param probabilities Vector of probabilities (must sum to 1)

PRIOR(PARAMDIST, DISTFACTORY) continuous form. @param paramDist Distribution of the scalar parameter theta @param distFactory Handle theta -> Distribution

@return self Prior instance

Property Summary

distFactory: Handle theta -> Distribution (continuous form)

distributions: Cell array of alternative distributions (discrete form)

kind: ‘discrete’ or ‘continuous’

paramDist: Distribution over the scalar parameter (continuous form)

probabilities: Vector of prior probabilities (sum to 1) (discrete form)

Method Summary

assertDiscrete(caller)

ASSERTDISCRETE(CALLER) Reject enumeration of a continuous prior.

getNumAlternatives returns NaN for a continuous prior, and every comparison against NaN is false, so a bounds check alone would pass and the caller would fault on an empty array instead.

discretize(n, method)

[DISTS, WEIGHTS] = DISCRETIZE(N, METHOD) Reduce the prior to N weighted alternatives.

The method is honoured for both forms of prior:

‘quadrature’ For a discrete prior, the alternatives and their: probabilities unchanged, and N is ignored: the set is already exact. For a continuous prior, stratified quantile midpoints with weights 1/N. Each node is the conditional median of an equal-mass stratum, so the rule integrates the parameter density in probability space and needs only evalCDF, which every Distribution provides.
‘montecarlo’ N i.i.d. draws, weights 1/N. For a discrete prior: the draws are of the alternative index against its probabilities; returning the alternatives unweighted here would silently drop the prior.

@param n Number of alternatives (ignored by a discrete quadrature) @param method ‘quadrature’ (default) or ‘montecarlo’ @return dists Cell array of Distribution objects @return weights Row vector of weights summing to 1

evalCDF(t)

FT = EVALCDF(T) Evaluate mixture CDF at t

F(t) = sum_i p_i * F_i(t)

evalLST(s)

L = EVALLST(S) Evaluate mixture Laplace-Stieltjes transform

L(s) = sum_i p_i * L_i(s)

static fromSample(k, s, distFactory)

SELF = FROMSAMPLE(K, S, DISTFACTORY) Continuous prior for a rate estimated from lifetime data.

Given K i.i.d. observations of an exponential random variable summing to S, the Jeffreys improper prior f(lambda) = s/lambda yields the posterior density of the rate

f(lambda|s) = lambda^(k-1) s^k exp(-lambda s) / (k-1)!

which is an Erlang density with K phases and phase rate S. See Trivedi and Bobbio (2017), Eq. (3.71). The posterior has mean K/S, i.e. the maximum-likelihood rate estimate, and variance K/S^2, so it concentrates on the estimate as K grows.

@param k Number of observations (positive integer) @param s Sum of the observed lifetimes (positive) @param distFactory Handle theta -> Distribution, default @(lambda) Exp(lambda) @return self Prior instance in continuous form

getAlternative(idx): DIST = GETALTERNATIVE(IDX) Return the distribution at index idx

getMean()

MEAN = GETMEAN() Get prior-weighted mean (expected mean over alternatives)

E[X] = sum_i p_i * E[X_i]

getNumAlternatives(): N = GETNUMALTERNATIVES() Return number of alternative distributions. A continuous prior has no alternatives until discretize() is called, so this returns NaN to force callers to discretize.

getProbability(idx): P = GETPROBABILITY(IDX) Return the probability of alternative idx

getSCV()

SCV = GETSCV() Get prior-weighted SCV using law of total variance

Var(X) = E[Var(X|D)] + Var(E[X|D]) SCV = Var(X) / E[X]^2

getSkewness(): SKEW = GETSKEWNESS() Get prior-weighted skewness (approximation using mixture formula)

isContinuousPrior(): BOOL = ISCONTINUOUSPRIOR() Return true if the prior is specified by a parameter density. Note: the base class isContinuous() refers to the support of the distribution itself, not to the form of the prior.

isPrior(): BOOL = ISPRIOR() Return true (used for detection by UQ solver)

static isPriorDistribution(dist)

BOOL = ISPRIORDISTRIBUTION(DIST) Check if a distribution is a Prior

@param dist Distribution object to check @return bool True if dist is a Prior

static quantile(dist, p)

X = QUANTILE(DIST, P) Numerical inverse CDF by bisection.

Uses only evalCDF, so it applies to any Distribution. Bracketing starts from the mean and doubles outward, which terminates for any distribution with finite mean.

@param dist Distribution object @param p Probability level in (0,1) @return x Value with F(x) = p

sample(n)

X = SAMPLE(N) Sample from prior (mixture sampling)

Samples are drawn from the mixture distribution where each sample comes from one of the alternatives selected according to the prior probabilities.

class Poisson

Poisson Discrete distribution for counting rare events

Poisson represents the Poisson distribution with rate parameter lambda. This distribution models the number of events occurring in a fixed interval when events occur independently at a constant average rate. It is commonly used for modeling arrival counts, defects, and other rare event phenomena.

@brief Poisson distribution for modeling discrete event counts

Key characteristics: - Single parameter: rate lambda (average number of events) - Mean = Variance = lambda - SCV = 1/lambda (decreases as lambda increases) - Support: {0, 1, 2, 3, …} - Limit of binomial distribution as n→∞, p→0, np→lambda

The Poisson distribution is used for: - Counting arrivals in fixed time intervals - Modeling defects or failures - Call center arrival modeling - Network packet arrival counts - Population dynamics and birth processes

Example: @code arrival_count = Poisson(3.5); % Average 3.5 arrivals per interval samples = arrival_count.sample(1000); % Generate counts @endcode

Constructor Summary

Poisson(lambda)

POISSON Create a Poisson distribution instance

@brief Creates a Poisson distribution with specified rate parameter @param lambda Rate parameter (average number of events, must be positive) @return self Poisson distribution instance

Method Summary

evalCDF(k): FT = EVALCDF(SELF,K) Evaluate the cumulative distribution function at K

evalLST(s): L = EVALST(S) Evaluate the Laplace-Stieltjes transform of the distribution function at s For Poisson(λ), LST(s) = exp(λ(e^(-s) - 1))

evalPMF(k): P = EVALPMF(K) Evaluate the probability mass function at K

getMean()

getProcess(): PROC = GETPROCESS()

getRate()

getSCV()

getSkewness(): Poisson skewness is lambda^(-1/2): E[(X-mu)^3] = lambda and sigma^3 = lambda^(3/2). The former sqrt(lambda) inverted it.

getVar()

sample(n): X = SAMPLE(N) Get n samples from the distribution

class Pareto

The Pareto statistical distribution

Constructor Summary

Pareto(shape, scale): SELF = PARETO(SHAPE, SCALE)

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T)

evalLST(s)

L = EVALST(S) Evaluate the Laplace-Stieltjes transform of the distribution function at s

A*(s) = E[e^{-sX}] = int_k^Inf e^{-sx} alpha k^alpha x^{-(alpha+1)} dx. Substituting x = k/u maps the infinite tail onto a unit interval and cancels the scale exactly:

A*(s) = alpha * int_0^1 u^(alpha-1) exp(-s*k/u) du

This is the same transform as the closed form of Nadarajah & Kotz, A*(s) = alpha*(s*k)^alpha*Gamma(-alpha, s*k) = alpha*E_{alpha+1}(s*k) (Queueing Syst (2006) 54:243-244, DOI 10.1007/s11134-006-0299-1), but in a form that stays accurate as s -> 0, where the incomplete-gamma product underflows to 0/Inf. Here s = 0 gives alpha*int_0^1 u^(alpha-1) du = 1 exactly, and the integrand is bounded and C^Inf on a FINITE interval for alpha >= 2 (the shape floor the constructor enforces).

Accuracy: adaptive Gauss-Kronrod at RelTol 1e-12, i.e. ~1e-12 relative, verified against mpmath to 1e-15. The previous implementation was a 1000-point right-endpoint rectangle sum truncated at k*1000^(1/alpha); it lost the mass beyond the truncation point and biased the transform low by ~3.1% at alpha=2.0078 (it returned A*(0)=0.96914, not 1).

static fitMeanAndSCV(MEAN, SCV): PA = FITMEANANDSCV(MEAN, SCV)

getMean(): EX = GETMEAN()

getProcess(): PROC = GETPROCESS()

getSCV(): SCV = GETSCV()

sample(n): X = SAMPLE(N)

class PH

Bases: Markovian

Abstract phase-type distribution with Markovian structure

Defined by initial probabilities and transient subgenerator matrix.

Constructor Summary

PH(alpha, T)

PH Create a phase-type distribution instance

@brief Creates a phase-type distribution with initial probabilities and subgenerator @param alpha Initial probability vector (must sum to ≤ 1) @param T Transient subgenerator matrix (must be nonsingular) @return self PH distribution instance

Method Summary

static fit(MEAN, SCV, SKEW): EX = FIT(MEAN, SCV, SKEW)

static fitCentral(MEAN, VAR, SKEW): EX = FITCENTRAL(MEAN, VAR, SKEW)

static fitMeanAndSCV(MEAN, SCV): EX = FITMEANANDSCV(MEAN, SCV)

static fitRawMoments(m1, m2, m3): Fit the distribution from first three moments

getInitProb(): ALPHA = GETINITPROB()

getSubgenerator(): T = GETSUBGENERATOR()

sample(n): X = SAMPLE(N)

setMean(MEAN): UPDATEMEAN(SELF,MEAN)

update(varargin): UPDATE(SELF,VARARGIN)

class NHPP

NHPP Non-homogeneous Poisson process with a piecewise-constant intensity.

The intensity is a step function of the wall clock: segment i covers [breakpoints(i), breakpoints(i+1)) and carries rate rates(i), so breakpoints has one more entry than rates. With cyclic=true the schedule repeats, giving a cyclic Poisson process.

Two horizon conventions:

cyclicthe schedule repeats with period: T = breakpoints(end) - breakpoints(1); the active segment at time t follows from mod(t - breakpoints(1), T).
non-cyclicthe intensity is zero outside: [breakpoints(1), breakpoints(end)), so the process emits nothing once the schedule is exhausted. A non-cyclic NHPP is therefore a transient construct: run to steady state it converges to the empty system, so callers should use a time span within the horizon.

This is NOT a renewal process. Successive intervals are dependent, because the position within the schedule carries over from one event to the next. Accordingly the scalar summaries that presuppose an i.i.d. interval distribution – getSCV, getSkewness, evalCDF, evalLST – are undefined and return NaN rather than a representative exponential value, which would silently misreport the process as Poisson. The schedule is the parameterisation: read it with getRateSchedule. getMean is well defined and returns the arrival-stationary (Palm) mean interval 1/timeAverageRate.

The process representation stores:: process{1} : 1-by-(n+1) row vector of breakpoints process{2} : 1-by-n row vector of rates process{3} : logical, true if cyclic

Solver support. The LDES simulation engine honours the exact schedule in both steady state (cyclic only) and transient analysis. SolverFLD honours it in getTranAvg, by injecting the intensity as a time-varying rate multiplier on the closing ODE; SolverFLD.getAvg uses the time-average rate, which is the steady state of a cyclic schedule. Every other solver rejects a model using it via the standard unsupported-feature check.

Constructor Summary

NHPP(breakpoints, rates, cyclic): SELF = NHPP(BREAKPOINTS, RATES, CYCLIC)

Property Summary

breakpoints: 1-by-(n+1) segment boundaries, strictly increasing

cyclic: logical, whether the schedule repeats

process: {breakpoints, rates, cyclic}, for serialization

rates: 1-by-n non-negative rate on each segment

sampleClock: wall-clock position of the next sample (see sample)

Method Summary

evalCDF(t): F = EVALCDF(T) NaN; see getSCV.

evalLST(s): L = EVALLST(S) NaN; see getSCV.

getBreakpoints()

getMean(): MEAN = GETMEAN() Arrival-stationary (Palm) mean interval.

getNumSegments()

getPeriod(): T = GETPERIOD() Horizon length, which is the period when cyclic.

getProcess()

getRateAt(t): R = GETRATEAT(T) Rate in force at T; zero past a non-cyclic horizon.

getRateSchedule(): SCHED = GETRATESCHEDULE() The parameterisation of the process; the scalar interval summaries are not. Model compilation recognises a schedule-bearing process by this method rather than by class name.

getRates()

getSCV(): SCV = GETSCV() NaN: an NHPP is not a renewal process, so there is no i.i.d. interval distribution for an SCV to summarise. Returning a representative value here would report a time-varying process as an exponential one to every consumer of sn.scv.

getSkewness(): SKEW = GETSKEWNESS() NaN; see getSCV.

getTimeAverageRate(): RATE = GETTIMEAVERAGERATE() sum(rates.*widths)/sum(widths). For a cyclic schedule this is the long-run arrival rate; for a non-cyclic one it averages over the active horizon only, the intensity being zero afterwards.

isCyclic()

nextInterval(from, residual)

X = NEXTINTERVAL(FROM, RESIDUAL) Solves int_{from}^{from+x} lambda(u) du = RESIDUAL for X, walking the schedule forward and consuming the budget segment by segment. Returns 0 when a non-cyclic horizon is exhausted first, which callers read as “no further event”.

Exact for an NHPP: conditional on no event since the last one, the residual is governed by the intensity from the current instant onward, so a holding time drawn under a rate that has since changed is not a sample from this process.

resetSampleClock(): RESETSAMPLECLOCK() Restart the sample path at the schedule start.

sample(n): X = SAMPLE(N) Draws N successive interarrival times along ONE sample path. The intensity depends on absolute time, so this advances an internal clock across calls: consecutive samples form a realisation of the process starting at breakpoints(1), not independent draws from a marginal. Use resetSampleClock to restart. A non-cyclic schedule that runs out returns 0 for every remaining sample, the intensity there being zero.

class MMPP2

2-phase Markov-Modulated Poisson Process - MMPP(2)

Constructor Summary

MMPP2(lambda0, lambda1, sigma0, sigma1): SELF = MMPP2(LAMBDA0,LAMBDA1,SIGMA0,SIGMA1)

Method Summary

D(i, wantSparse): Di = D(i)

evalACFT(lags, timescale)

ACF = EVALACFT(self, lags)

Evaluate the autocorrelation in counts at timescale t

evalMeanT(t): MEANT = EVALMEANT(SELF,T)

evalVarT(t): VART = EVALVART(SELF,T)

static fitCentralAndACFDecay(mean, var, skew, gamma2)

static fitCentralAndACFLag1(mean, var, skew, rho1)

static fitCentralAndIDC(mean, var, skew, idc)

static fitRawMomentsAndACFDecay(m1, m2, m3, gamma2)

static fitRawMomentsAndACFLag1(m1, m2, m3, rho1)

static fitRawMomentsAndIDC(m1, m2, m3, idc)

getACF(lags): ACF = GETACF(self, lags)

getACFDecay()

[gamma2, gamma] = GETACFDECAY(self)

gamma2: asymptotic decay rate of acf gamma: interpolated decay rate of acf

getIDC(): IDC = GETIDC() % INDEX OF DISPERSION FOR COUNTS

getIDI(): IDI = GETIDI() % INDEX OF DISPERSION FOR INTERVALS

getMean(): MEAN = GETMEAN()

getNumberOfPhases(~): N = GETNUMBEROFPHASES()

getRate()

getSCV(): SCV = GETSCV()

getSkewness()

isImmediate(): BOOL = ISIMMEDIATE()

static rand()

MMPP2 = RAND()

Generate random MAP using uniform random numbers

sample(n): X = SAMPLE(N)

toMMDP()

TOMMDP Convert MMPP2 to MMDP2 (deterministic representation)

Converts this Markov-Modulated Poisson Process to a Markov-Modulated Deterministic Process suitable for fluid queue analysis.

@return mmdp MMDP2 representation of this MMPP2

class ME

Bases: Markovian

Matrix Exponential (ME) distribution

ME distributions are characterized by an initial vector alpha and a matrix parameter A. They generalize Phase-Type (PH) distributions by allowing alpha to have entries outside [0,1] and A to have arbitrary structure (not necessarily a valid sub-generator).

Constructor Summary

ME(alpha, A, checkDensity)

ME Create a Matrix Exponential distribution instance

@brief Creates an ME distribution with the given initial vector and matrix parameter @param alpha Initial vector (may have negative entries or sum != 1) @param A Matrix parameter (must have all eigenvalues with negative real parts) @param checkDensity Scan the density for a negative value (default true).

Set to false only by subclasses whose representation is a density by construction, such as CME, where the scan cannot fire and costs O(1e5) propagations of a large matrix.

@return self ME distribution instance

Property Summary

A: Matrix parameter

alpha: Initial vector

Method Summary

evalCDF(t)

FT = EVALCDF(SELF, T) Evaluate the cumulative distribution function at t

For ME distribution: CDF(t) = 1 - alpha * exp(A*t) * e

evalLST(s): L = EVALST(S) Evaluate the Laplace-Stieltjes transform at s

evalPDF(t)

FT = EVALPDF(SELF, T) Evaluate the probability density function at t

For ME distribution: PDF(t) = -alpha * exp(A*t) * A * e

static fitMoments(moms)

ME = FITMOMENTS(MOMS) Create ME distribution by fitting the given moments Uses BuTools MEFromMoments algorithm

@param moms Array of moments (requires 2*M-1 moments for order M) @return me ME distribution matching the given moments

static fromErlang(k, rate)

ME = FROMERLANG(K, RATE) Create ME distribution from Erlang distribution Convenience method showing that Erlang is a special case of ME

@param k Number of phases @param rate Rate parameter for each phase @return me ME distribution equivalent to Erlang(k, rate)

static fromExp(rate)

ME = FROMEXP(RATE) Create ME distribution from exponential distribution Convenience method showing that Exp is a special case of ME

@param rate Rate parameter (lambda) @return me ME distribution equivalent to Exp(rate)

static fromHyperExp(p, rates)

ME = FROMHYPEREXP(P, RATES) Create ME distribution from HyperExponential distribution Convenience method showing that HyperExp is a special case of ME

@param p Array of probabilities for each branch @param rates Array of rates for each branch @return me ME distribution equivalent to HyperExp(p, rates)

getMean()

MEAN_VAL = GETMEAN() Get mean of the ME distribution, m1 = -alpha*inv(A)*e

The definition is used rather than map_mean, which obtains the rate from the stationary vector of D0+D1. That vector is a probabilistic object of a Markovian process, and solving for it on an ME degrades with the oscillation of A: a CME of order 101 came out with a relative error of 5.6e-8 in the mean and 2.6e-4 in the SCV, where the definition is exact to 1e-13. Native Python ME.getMean and jline.lang.processes.ME do the same.

getNumberOfPhases(): PHASES = GETNUMBEROFPHASES() Get number of phases in the ME representation

getProcess(): PROC = GETPROCESS() Get process representation {D0, D1}

getSCV(): SCV = GETSCV() Get squared coefficient of variation, var/mean^2

getVar(): VAR_VAL = GETVAR() Get variance of the ME distribution, m2 - m1^2 with m2 = 2*alpha*inv(A)^2*e

sample(n): X = SAMPLE(N) Get n samples from the distribution using inverse CDF interpolation

static scanNegativeDensity(alpha, A)

[ISNEG, FMIN, TMIN] = SCANNEGATIVEDENSITY(ALPHA, A) Search the density f(t) = -alpha*expm(A*t)*A*e for a negative value.

A negative value found here is a witness: it proves that the representation is not a distribution. Finding none proves nothing, so the caller must not report the converse.

This replaces CheckMEPositiveDensity as the trigger for the construction-time warning. That routine searches for a Markovian monocyclic equivalent, which is a sufficient condition only, and its verdict depends on the representation rather than on the distribution: for alpha=[1,0,0], A=[-0.5 0 0; 0 -1 w; 0 -w -1] the distribution is Exp(0.5) for every w, yet the search fails once w >= 2*pi. It also costs of the order of a second per call at search order 1000, which is far too slow for a constructor.

The horizon covers all but scanTail of the mass, using the dominant (least negative) eigenvalue of A; the sampling rate resolves the fastest oscillation present, taken from the largest imaginary part. The constants match native Python line_solver.distributions.markovian and jline.lang.processes.ME.

class MAP

Markovian Arrival Process for correlated arrival modeling

Models arrival streams with correlation and burstiness via D0 and D1 matrices.

Constructor Summary

MAP(D0, D1)

MAP Create a Markovian Arrival Process instance

@brief Creates a MAP with specified D0 and D1 matrices @param D0 Generator matrix for transitions without arrivals @param D1 Rate matrix for transitions with arrivals @return self MAP instance with specified matrices

Method Summary

D(i, wantSparse): Di = D(i)

evalACFT(lags, timescale)

ACF = EVALACFT(self, lags)

Evaluate the autocorrelation in counts at timescale t

evalMeanT(t): MEANT = EVALMEANT(SELF,T)

evalVarT(t): VART = EVALVART(SELF,T)

getACF(lags): ACF = GETACF(self, lags)

getACFDecay()

[gamma2, gamma] = GETACFDECAY(self)

gamma2: asymptotic decay rate of acf gamma: interpolated decay rate of acf

getIDC(t): index of dispersion for counts IDC = GETIDC() % ASYMPTOTIC INDEX OF DISPERSION

getNumberOfPhases(): The method call needs its receiver: without it this reads as a call to a function named D and errors out, which made every model holding a MAP fail at getStruct time (refreshPetriNetNodes.m:50 asks a Markovian firing distribution for its phase count).

getRate(): MAP = GETRATE()

isImmediate()

static rand(order)

MAP = RAND(ORDER)

Generate random MAP using uniform random numbers

static randn(order, mu, sigma)

MAP = RANDN(ORDER, MU, SIGMA)

Generate random MAP using specified Gaussian parameter and taking the absolute value of the resulting values

sample(n): X = SAMPLE(N)

setMean(MEAN): UPDATEMEAN(SELF,MEAN) Update parameters to match the given mean

toMMDP()

TOMMDP Convert MAP to MMDP (deterministic representation)

Converts this Markovian Arrival Process to a Markov-Modulated Deterministic Process suitable for fluid queue analysis.

@return mmdp MMDP representation of this MAP

toTimeReversed()

class Geometric

A Geometric probability distribution

The distribution of the number of Bernoulli trials needed to get one success.

Constructor Summary

Geometric(p): SELF = GEOMETRIC(P)

Method Summary

evalCDF(k): FT = EVALCDF(SELF,K)

evalLST(s): L = EVALST(S) Evaluate the Laplace-Stieltjes transform of the distribution function at s For Geometric(p), LST(s) = p*e^(-s) / (1 - (1-p)*e^(-s))

evalPMF(k): PR = EVALPMF(K)

getMean(): EX = GETMEAN()

getProcess(): PROC = GETPROCESS()

getSCV(): SCV = GETSCV()

sample(n): X = SAMPLE(N)

class Cox2

Bases: Markovian

Static class to fit two-phase coxian statistical distribution

Method Summary

static fit(MEAN, SCV, SKEW): CX = FIT(MEAN,SCV,SKEW)

static fitCentral(MEAN, VAR, SKEW): CX = FITCENTRAL(MEAN,VAR,SKEW)

static fitMean(MEAN): CX = FITMEAN(MEAN)

static fitMeanAndSCV(MEAN, SCV): CX = FITMEANANDSCV(MEAN, SCV)

class CME

Bases: ME

Concentrated Matrix Exponential (CME) distribution

A CME is the matrix-exponential distribution of odd order 2*n+1 whose squared coefficient of variation is (numerically) minimal for that order, from the tables of Horvath, Horvath and Telek. Its SCV decays as O(1/n^2) and therefore goes far below the Erlang bound 1/order attainable by a phase-type distribution of the same order: order 101 gives SCV 3.9e-4, where Erlang-101 gives 9.9e-3.

The density of the unit-mean CME with n harmonic terms is

f(x) = mu1*exp(-mu1*x)*(c + sum_k [a_k*cos(k*w*mu1*x) + b_k*sin(k*w*mu1*x)])

with w = omega, which is exactly alpha*expm(A*x)*(-A*e) for the block-diagonal A = blkdiag(-mu1, mu1*[-1 -k*w; k*w -1], k=1..n). The parameters a, b, c, omega and mu1 are read from the same iltcme.json table used by the CME inverse Laplace transform (matlab_ilt).

Constructor Summary

CME(mean, order)

SELF = CME(MEAN, ORDER) Create a concentrated matrix-exponential distribution

@param mean Mean of the distribution (positive) @param order Number of phases, an odd integer 2*n+1 with n in the table @return self CME distribution instance

Property Summary

cmeMean: Mean of the distribution

cmeOrder: Number of phases, an odd integer 2*n+1

Method Summary

static fitMeanAndSCV(mean, scv)

CME = FITMEANANDSCV(MEAN, SCV) Create the lowest-order CME with the given mean and SCV at most scv

@param mean Target mean @param scv Target squared coefficient of variation, an upper bound.

The lowest tabulated order whose minimal SCV does not exceed it is selected, so the result is at least as concentrated as requested.

@return cme CME of that order rescaled to the requested mean

static getMinSCV(order): SCV = GETMINSCV(ORDER) Get the tabulated minimal SCV attained by a CME of the given order

getOrder(): O = GETORDER() Get the CME order, i.e. the number of phases

static getSupportedOrders(): ORDERS = GETSUPPORTEDORDERS() Get the sorted list of CME orders (phase counts) in the table

static representation(order)

[ALPHA, A, SCV] = REPRESENTATION(ORDER) Build the unit-mean (alpha, A) matrix-exponential form of a CME

A = blkdiag(-mu1, mu1*[-1 -k*w; k*w -1], k=1..n) reproduces the exponential envelope in its first phase and the k-th harmonic in its k-th 2x2 rotation block, since expm(mu1*x*[-1 -kw; kw -1]) is exp(-mu1*x) times the rotation by k*w*mu1*x. The entries of alpha follow by matching -alpha*expm(A*x)*A*e term by term: with wk = k*w and d = 2*(1+wk^2),

alpha(1) = c alpha(2k) = ((1+wk)*a_k - (1-wk)*b_k)/d alpha(2k+1) = ((1-wk)*a_k + (1+wk)*b_k)/d

The result has unit mean and sums to one, as any (alpha, A) whose density integrates to one must.

static table()

PARAMS = TABLE() Load and cache the CME parameter table from iltcme.json

The same table backs matlab_ilt, which caches it in the global cmeParams, so the two share one decode per session.

static tableEntry(order)

ENTRY = TABLEENTRY(ORDER) Select the CME table entry realizing the given number of phases

The table is keyed by the number of harmonic terms n, so an order of 2*n+1 phases maps to the entries with that n. Several entries can share an n (the ‘full’ and ‘approx’ optimizations), and the most concentrated one is taken, matching the selection rule of the CME inverse Laplace transform.

class BMAP

Bases: MarkedMAP

Batch Markovian Arrival Process (BMAP)

BMAP is a point process where arrivals occur in batches. Uses the standard BMAP representation: - D0: infinitesimal generator for transitions without arrivals - D1: rate matrix for transitions generating 1 arrival - D2: rate matrix for transitions generating 2 arrivals - … - Dk: rate matrix for transitions generating k arrivals

BMAP extends MarkedMAP where each “mark” k represents a batch size k.

Constructor Summary

BMAP(D)

SELF = BMAP(D)

D is a cell array {D0, D1, D2, …, Dk} where: - D0: transitions without arrivals - Dk: transitions generating k arrivals (for k >= 1)

The number of marking types K equals the maximum batch size (i.e., K = length(D) - 2 when D also includes D1_total,

or K = length(D) - 1 otherwise)

Method Summary

display(): DISPLAY(SELF)

static fromMAPWithBatchPMF(D0, D1, batchSizes, pmf)

BMAP = FROMMAPWITHBATCHPMF(D0, D1, BATCHSIZES, PMF)

Create BMAP from a base MAP and batch size distribution

Inputs:: D0: base MAP’s D0 matrix (transitions without batch arrivals) D1: base MAP’s D1 matrix (inter-batch arrival transitions) batchSizes: array of batch sizes (e.g., [1, 2, 4, 8]) pmf: probability mass function for batch sizes (must sum to 1)
Output:: bmap: BMAP constructed from the base MAP and batch distribution

getBatchRates()

RATES = GETBATCHRATES()

Returns array where rates(k) is the rate of batch size k arrivals

getInterBatchMAP()

MAP = GETINTERBATCHMAP()

Returns the underlying MAP for inter-batch arrivals

getMaxBatchSize()

K = GETMAXBATCHSIZE()

Returns the maximum batch size k where Dk is defined

getMeanBatchSize()

MEAN_BS = GETMEANBATCHSIZE()

Computes the mean batch size as a weighted average: E[batch size] = sum(k * rate_k) / sum(rate_k)

getName()

NAME = GETNAME()

Returns the name of this distribution type

static rand(order, maxBatchSize)

BMAP = RAND(ORDER, MAXBATCHSIZE)

Generate random BMAP using uniform random numbers

Inputs:: order: number of phases (default: 2) maxBatchSize: maximum batch size (default: 3)

sample(n)

[X, B] = SAMPLE(N)

Sample n batches from the BMAP :returns: X – inter-arrival times (between batches)

B: batch sizes for each arrival

sampleInterBatch(n)

X = SAMPLEINTERBATCH(N)

Sample n inter-batch arrival times (ignoring batch sizes)

validateGenerator()

class APH

Bases: Markovian

APH Acyclic Phase-type distribution with tree-like structure

APH represents acyclic phase-type distributions where the underlying Markov chain has no cycles (tree-like structure). This subclass of phase-type distributions provides computational advantages while maintaining flexibility for modeling various service time distributions.

@brief Acyclic phase-type distribution with tree-like Markov structure

Key characteristics: - Subclass of general phase-type distributions - Acyclic structure (no cycles in transition graph) - Tree-like Markov chain topology - Computational advantages over general PH distributions - Maintains distribution approximation capabilities

Acyclic structure benefits: - More efficient algorithms for many operations - Simpler parameter fitting procedures - Easier interpretation and visualization - Reduced computational complexity in analyses

APH distributions are used for: - Service time modeling with reduced complexity - Fitting distributions with controlled structure - Matrix-analytic methods requiring efficiency - Systems where acyclic behavior is natural - Performance models needing fast computation

Example: @code alpha = [1, 0]; T = [-2, 1; 0, -3]; % Simple 2-phase acyclic aph_dist = APH(alpha, T); samples = aph_dist.sample(1000); @endcode

Constructor Summary

APH(alpha, T)

APH Create an acyclic phase-type distribution instance

@brief Creates an APH distribution with acyclic structure @param alpha Initial probability vector or Java APH object @param T Transient subgenerator matrix (acyclic structure) @return self APH distribution instance

Method Summary

evalCDF(varargin)

static fit(MEAN, SCV, SKEW): EX = FIT(MEAN, SCV, SKEW)

static fitCentral(MEAN, VAR, SKEW): EX = FITCENTRAL(MEAN, VAR, SKEW)

static fitMeanAndSCV(MEAN, SCV): EX = FITMEANANDSCV(MEAN, SCV)

static fitRawMoments(m1, m2, m3): Fit the distribution from first three moments

getInitProb(): ALPHA = GETINITPROB()

getSubgenerator(): T = GETSUBGENERATOR()

static rand(order)

sample(n): X = SAMPLE(N)

setMean(MEAN): UPDATEMEAN(SELF,MEAN)

update(varargin): UPDATE(SELF,VARARGIN)

updateFromRawMoments(varargin): UPDATE(SELF,VARARGIN)

class MarkedMAP

Markov Modulated Arrival Process

Constructor Summary

MarkedMAP(D, K)

SELF = MarkedMAP(D,K)

LINE uses the M3A representation format D={D0, D1, D11, D12, D13, …, D1K} K is the number of marking types

Method Summary

D(i, j, wantSparse): Di = D(i)

evalACFT(lags, timescale)

ACF = EVALACFT(self, lags)

Evaluate the autocorrelation in counts at timescale t

evalMeanT(t): MEANT = EVALMEANT(SELF,T)

evalVarT(t): VART = EVALVART(SELF,T)

static fit(trace, markings, order): M3PP = FIT(TRACE, ORDER)

getACF(lags): ACF = GETACF(self, lags)

getACFDecay()

[gamma2, gamma] = GETACFDECAY(self)

gamma2: asymptotic decay rate of acf gamma: interpolated decay rate of acf

getIDC(t): index of dispersion for counts ID = GETIDC() % INDEX OF DISPERSION

getMarkedIDC(t): asymptotic index of dispersion for counts for each type ID = GETMARKEDIDC() % ASYMPTOTIC INDEX OF DISPERSION

getMarkedMeans(): MEAN = GETMARKEDMEANS()

getMean(): MEAN = GETMEAN()

getMu(): MU = GETMU() Aggregate departure rate from each state

getNumberOfPhases(): N = GETNUMBEROFMAPASES()

getNumberOfTypes(): K = GETNUMBEROFTYPES() Number of marking types

getPhi(): MAPI = GETMAPI() Return the exit vector of the underlying MAP

getRate(): LAMBDA = GETRATE()

getSCV(): SCV = GETSCV()

static rand(order, nclasses)

MarkedMAP = RAND(ORDER,NCLASSES)

Generate random MarkedMAP using uniform random numbers

sample(n): [X,C] = SAMPLE(N)

toMAP()

toMAPs(types)

toTimeReversed()

class Marked

An abstract class for Markov-modulated marked processes

Method Summary

getEmbedded()

P = GETEMBEDDED()

Get DTMC embedded at event arrival times

getEmbeddedProb()

PIE = GETEMBEDDEDRPOB()

Solve DTMC embedded embedded at event arrival times

sample(n): X = SAMPLE(N)

class Zipf

Zipf Power-law distribution for modeling popularity and ranking

Zipf represents the Zipf distribution commonly used for modeling popularity distributions, where items are ranked by frequency. The probability of an item is inversely proportional to its rank raised to the shape parameter. This distribution is essential for cache modeling and content popularity.

@brief Zipf distribution for popularity and ranking-based phenomena

Key characteristics: - Two parameters: shape (s) and number of items (n) - Probability ∝ 1/rank^s for each rank - Heavy-tailed distribution (popular items dominate) - Mean = H(s-1,n)/H(s,n) where H is generalized harmonic number - Support: {1, 2, …, n}

Shape parameter behavior: - s > 1: Heavy tail, few items very popular - s = 1: Classic Zipf law - s < 1: Less skewed, more uniform popularity

The Zipf distribution is used for: - Web cache popularity modeling - Content access patterns - Word frequency in natural language - City population distributions - Internet traffic modeling

Example: @code web_popularity = Zipf(1.2, 1000); % 1000 items, shape=1.2 access_pattern = web_popularity.sample(10000); @endcode

Constructor Summary

Zipf(s, n)

ZIPF Create a Zipf distribution instance

@brief Creates a Zipf distribution with shape parameter and item count @param s Shape parameter (s > 0, controls skewness) @param n Number of items (positive integer) @return self Zipf distribution instance with specified parameters

Method Summary

evalCDF(k): FT = EVALCDF(SELF,K)

evalLST(s): L = EVALST(S) Evaluate the Laplace-Stieltjes transform of the distribution function at s For Zipf distribution, LST(s) = Σ(k=1 to n) P(X=k) * e^(-s*k)

evalPMF(k): P = EVALPMF(K)

static genHarmonic(s, n): HNM = GENHARMONIC(S,N)

getMean(): EX = GETMEAN()

getProcess(): PROC = GETPROCESS()

getSCV(): SCV = GETSCV()

sample(n): X = SAMPLE(N) to be checked

class Weibull

Weibull Distribution for reliability and failure time modeling

Weibull represents the Weibull distribution with shape and scale parameters. This distribution is widely used in reliability engineering and survival analysis for modeling failure times, life distributions, and extreme value phenomena. It can model increasing, decreasing, or constant hazard rates.

@brief Weibull distribution for reliability and failure time analysis

Key characteristics: - Two parameters: scale (alpha) and shape (r) - Mean = alpha * Γ(1 + 1/r) where Γ is the gamma function - Flexible hazard rate behavior based on shape parameter - Support: (0, ∞) - Includes exponential (r=1) as special case

Hazard rate behavior: - r < 1: Decreasing failure rate (infant mortality) - r = 1: Constant failure rate (exponential, random failures) - r > 1: Increasing failure rate (wear-out failures)

The Weibull distribution is used for: - Reliability and survival analysis - Failure time modeling - Lifetime data analysis - Wind speed modeling - Service time distributions with varying hazard rates

Example: @code failure_dist = Weibull(2.0, 1000); % Scale=1000, Shape=2 (wear-out) reliability_time = failure_dist.sample(100); @endcode

Constructor Summary

Weibull(shape, scale)

WEIBULL Create a Weibull distribution instance

@brief Creates a Weibull distribution with specified shape and scale @param shape Shape parameter r (must be positive) @param scale Scale parameter alpha (must be positive) @return self Weibull distribution instance

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T)

evalLST(s): L = EVALST(S) Evaluate the Laplace-Stieltjes transform of the distribution function at s The LST of Weibull distribution doesn’t have a simple closed form. We compute it numerically using the definition: E[e^(-sX)] = ∫₀^∞ e^(-sx) f(x) dx where f(x) is the Weibull PDF

static fitMeanAndSCV(MEAN, SCV): PA = FITMEANANDSCV(MEAN, SCV)

getMean(): EX = GETMEAN()

getProcess(): PROC = GETPROCESS()

getSCV(): SCV = GETSCV()

sample(n): X = SAMPLE(N)

class Process

Bases: Copyable

An abstract class for stochastic processes

Property Summary

name

params

class Normal

Normal Normal (Gaussian) distribution

Normal represents a normal (Gaussian) distribution with mean mu and standard deviation sigma. This distribution is widely used in modeling continuous phenomena due to the central limit theorem.

@brief Normal distribution with specified mean and standard deviation

Key characteristics: - Two parameters: mean (mu) and standard deviation (sigma) - Variance = sigma^2 - SCV = sigma^2 / mu^2 - Support: (-Inf, Inf) - Symmetric with zero skewness

Note: For queueing applications, a truncated or shifted version may be needed since normal distributions can take negative values.

Example: @code dist = Normal(5.0, 1.0); % Mean=5, StdDev=1 samples = dist.sample(1000); @endcode

Constructor Summary

Normal(mu, sigma)

NORMAL Create a Normal distribution instance

@brief Creates a Normal distribution with specified mean and std dev @param mu Mean of the distribution @param sigma Standard deviation of the distribution (must be positive) @return self Normal distribution instance

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T)

evalLST(s): L = EVALST(S)

static fitMean(MEAN): NM = FITMEAN(MEAN)

static fitMeanAndSCV(MEAN, SCV): NM = FITMEANANDSCV(MEAN, SCV)

static fitMeanAndStd(MEAN, STD): NM = FITMEANANDSTD(MEAN, STD)

static fitMeanAndVar(MEAN, VAR): NM = FITMEANANDVAR(MEAN, VAR)

getMean(): EX = GETMEAN()

getProcess(): PROC = GETPROCESS()

getSCV(): SCV = GETSCV()

getSkewness(): SKEW = GETSKEWNESS()

getStd(): STD = GETSTD()

getVar(): VAR = GETVAR()

sample(n): X = SAMPLE(N)

class MultivariateNormal

MultivariateNormal Multivariate Normal (Gaussian) Distribution

Represents a d-dimensional normal distribution with mean vector mu and covariance matrix Sigma. The distribution can be used standalone or within a Prior for mixture models.

Constructor:: mvn = MultivariateNormal(mu, Sigma)

Parameters:

mu - d-dimensional mean vector (column vector)
Sigma - d x d positive definite covariance matrix

Examples

% Create 2D normal distribution mu = [1; 2]; Sigma = [1, 0.5; 0.5, 1]; mvn = MultivariateNormal(mu, Sigma);

% Generate samples samples = mvn.sample(100); % 100 x 2 matrix

% Evaluate PDF pdf_val = mvn.evalPDF([1; 2]);

% Extract marginal norm = mvn.getMarginalUniv(1);

Constructor Summary

MultivariateNormal(mu, Sigma): Validate inputs

Property Summary

dimension: Dimensionality

Method Summary

evalCDF(t): CDF is not well-defined for multivariate distributions

evalLST(s): LST is not well-defined for multivariate distributions

evalPDF(x)

Evaluate the multivariate normal PDF

Parameters:: x - d x 1 column vector or d x n matrix of points
Returns:: p - PDF value(s)

static fitMeanAndCovariance(mu, Sigma)

Create a multivariate normal from mean and covariance

Parameters:

mu - mean vector
Sigma - covariance matrix

Returns:

mvn - MultivariateNormal distribution

getCorrelation(): Get the correlation matrix

getCovariance(): Get the covariance matrix (d x d)

getDimension(): Get the dimensionality of the distribution

getMarginal(indices)

Extract a marginal distribution for a subset of dimensions

Parameters:: indices - vector of dimension indices to keep (1-based)
Returns:: mvn_marg - MultivariateNormal for the marginal

getMarginalUniv(index)

Extract a univariate marginal distribution

Parameters:: index - dimension index (1-based)
Returns:: norm - 1D MultivariateNormal distribution for that dimension

getMean(): Get the mean of the first component (for Prior compatibility)

getMeanVector(): Get the mean vector (d x 1)

getSCV(): Get squared coefficient of variation (not meaningful for multivariate)

getSkewness(): Get skewness (multivariate normal is symmetric)

getVar(): Get the variance of the first component

sample(n)

Generate n samples from the multivariate normal distribution

Returns n x d matrix of samples

toString(): Convert to string representation

class MarkovModulated

Bases: Markovian

An abstract class for Markov-modulated processes

Method Summary

getEmbedded()

P = GETEMBEDDED()

Get DTMC embedded at event arrival times

getEmbeddedProb()

PIE = GETEMBEDDEDRPOB()

Solve DTMC embedded embedded at event arrival times

sample(n): X = SAMPLE(N)

class MarkedMarkovProcess

Bases: MarkovProcess

A class for continuous time Markov chain where transitions are labeled

Constructor Summary

MarkedMarkovProcess(InfGen, eventFilt, evs, isFinite, stateSpace): SELF = MARKEDMARKOVPROCESS(InfGen, eventFilt, events, isFinite, stateSpace)

Property Summary

eventFilt

eventList

Method Summary

embeddedSolve(evset): PIE = EMBEDDEDSOLVE(EVENTSET)

static fromSampleSysAggr(sa)

toMAP(varargin): MAP = TOMAP(EVENT)

class MarkedMMPP

Marked Markov-Modulated Poisson Process Also referred to as a M3PP model

Constructor Summary

MarkedMMPP(D, K)

SELF = MarkedMMPP(D,K)

LINE uses the M3A representation format D={D0, D1, D11, D12, D13, …, D1K} K is the number of marking types

Method Summary

D(i, j, wantSparse): Di = D(i)

evalACFT(lags, timescale)

ACF = EVALACFT(self, lags)

Evaluate the autocorrelation in counts at timescale t

evalMeanT(t): MEANT = EVALMEANT(SELF,T)

evalVarT(t): VART = EVALVART(SELF,T)

static fit(trace, markings, order): M3PP = FIT(TRACE, ORDER)

getACF(lags): ACF = GETACF(self, lags)

getACFDecay()

[gamma2, gamma] = GETACFDECAY(self)

gamma2: asymptotic decay rate of acf gamma: interpolated decay rate of acf

getIDC(t): index of dispersion for counts ID = GETIDC() % INDEX OF DISPERSION

getMarkedIDC(t): asymptotic index of dispersion for counts for each type ID = GETMARKEDIDC() % ASYMPTOTIC INDEX OF DISPERSION

getMarkedMeans(): MEAN = GETMARKEDMEAN()

getMean(): MEAN = GETMEAN()

getMu(): MU = GETMU() Aggregate departure rate from each state

getNumberOfPhases(): N = GETNUMBEROFMAPASES()

getNumberOfTypes(): K = GETNUMBEROFTYPES() Number of marking types

getPhi(): MAPI = GETMAPI() Return the exit vector of the underlying MAP

getProcess(): MarkedMMPP = GETPROCESS()

getRate(): LAMBDA = GETRATE()

getSCV(): SCV = GETSCV()

static rand(order, nclasses)

M3PP = RAND(ORDER,NCLASSES)

Generate random MarkedMMPP using uniform random numbers

sample(n): [X,C] = SAMPLE(N)

toMAP()

toMAPs(types)

toTimeReversed()

class MMDP2

Bases: MMDP

2-state Markov-Modulated Deterministic Process

A specialized MMDP with exactly 2 phases, using a convenient parameterization analogous to MMPP2.

Parameterization:: r0, r1: Deterministic rates in states 0 and 1 sigma0: Transition rate from state 0 to state 1 sigma1: Transition rate from state 1 to state 0
The generator matrix is:: Q = [-sigma0, sigma0; sigma1, -sigma1]
The rate matrix is:: R = diag([r0, r1])

Constructor Summary

MMDP2(r0, r1, sigma0, sigma1)

MMDP2 Create a 2-state Markov-Modulated Deterministic Process

@brief Creates a 2-state MMDP with specified rates and transitions @param r0 Deterministic rate in state 0 @param r1 Deterministic rate in state 1 @param sigma0 Transition rate from state 0 to state 1 @param sigma1 Transition rate from state 1 to state 0 @return self MMDP2 instance

Method Summary

Q()

Q Return the generator matrix

@return Q_mat 2×2 generator matrix

R()

R Return the rate matrix (diagonal)

@return R_mat 2×2 diagonal rate matrix

getMeanRate()

GETMEANRATE Compute stationary mean rate (closed-form)

For a 2-state MMDP, the mean rate has the closed form:: E[r] = (r0*σ1 + r1*σ0) / (σ0 + σ1)

@return rate Stationary mean deterministic rate

getNumberOfPhases(~)

GETNUMBEROFPHASES Return the number of phases

@return n Always 2 for MMDP2

getSCV()

GETSCV Compute squared coefficient of variation (closed-form)

For a 2-state MMDP, the SCV has a closed form based on the variance of rates over the stationary distribution.

@return scv Squared coefficient of variation

r()

r Return the rate vector

@return r_vec 2-vector of rates [r0; r1]

class MMDP

Bases: Markovian

Markov-Modulated Deterministic Process for fluid queue modeling

Models fluid flow with deterministic rates modulated by a background Markov chain. Suitable for arrival and service processes in Markovian fluid queues analyzed by the mfq method of SolverFLD.

The (Q, R) parameterization follows BUTools conventions: - Q: Generator matrix of the modulating CTMC - R: Diagonal matrix of deterministic rates per state

MMDP is the deterministic analogue of MMPP: - MMPP: Poisson arrival rates modulated by a Markov chain - MMDP: Deterministic rates modulated by a Markov chain

Constructor Summary

MMDP(Q, R)

MMDP Create a Markov-Modulated Deterministic Process

@brief Creates an MMDP with specified generator Q and rate matrix R @param Q n×n generator matrix (row sums = 0) @param R n×n diagonal matrix of rates, OR n-vector @return self MMDP instance with specified matrices

Method Summary

Q()

Q Return the generator matrix

@return Q_mat n×n generator matrix of the modulating CTMC

R()

R Return the rate matrix (diagonal)

@return R_mat n×n diagonal matrix of deterministic rates

static fromMAP(map_obj)

FROMMAP Convert a MAP to MMDP (deterministic representation)

Converts a Markovian Arrival Process to a Markov-Modulated Deterministic Process by extracting the full generator and using row sums of D1 as the deterministic rates.

@param map_obj MAP object to convert @return mmdp MMDP representation of the MAP

static fromMMPP2(lambda0, lambda1, sigma0, sigma1)

FROMMMPP2 Create MMDP from MMPP2 parameters

Creates a 2-state MMDP using the same parameterization as MMPP2.

@param lambda0 Rate in state 0 @param lambda1 Rate in state 1 @param sigma0 Transition rate from state 0 to state 1 @param sigma1 Transition rate from state 1 to state 0 @return mmdp 2-state MMDP

getMean()

GETMEAN Return mean inter-arrival time (inverse of mean rate)

@return mean Mean inter-arrival time (1/rate), or Inf if rate=0

getMeanRate()

GETMEANRATE Compute stationary mean rate

Computes E[r] = π * r, where π is the stationary distribution of the modulating CTMC with generator Q.

@return rate Stationary mean deterministic rate

getNumberOfPhases()

GETNUMBEROFPHASES Return the number of phases

@return n Number of phases in the modulating CTMC

getRate()

GETRATE Alias for getMeanRate

@return rate Stationary mean deterministic rate

getSCV()

GETSCV Return squared coefficient of variation of rates

Computes Var[r]/E[r]^2 where expectation is over the stationary distribution of the modulating CTMC.

@return scv Squared coefficient of variation

isImmediate()

ISIMMEDIATE Check if the process has infinite rate

@return bool True if mean rate is effectively infinite

r()

r Return the rate vector (diagonal of R)

@return r_vec n-vector of deterministic rates per phase

class Lognormal

The Lognormal statistical distribution

Constructor Summary

Lognormal(mu, sigma): SELF = LOGNORMAL(MU, SIGMA)

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T)

evalLST(s): L = EVALST(S) Evaluate the Laplace-Stieltjes transform of the distribution function at s The LST of Lognormal distribution doesn’t have a simple closed form. We compute it numerically using the definition: E[e^(-sX)] = ∫₀^∞ e^(-sx) f(x) dx where f(x) is the lognormal PDF

static fitMeanAndSCV(MEAN, SCV): PA = FITMEANANDSCV(MEAN, SCV)

getMean(): EX = GETMEAN()

getProcess(): PROC = GETPROCESS()

getSCV(): SCV = GETSCV()

sample(n): X = SAMPLE(N)

class Immediate

Bases: Distribution

A distribution with probability mass entirely at zero

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T)

evalLST(s): L = EVALST(S)

static getInstance()

getMean(): EX = GETMEAN()

getMu(): MU = GETMU() global GlobalConstants.Immediate

Return total outgoing rate from each state

getPhi(): PHI = GETPHI()

getProcess()

getRate(): global GlobalConstants.Immediate

getSCV(): SCV = GETSCV()

isDisabled(): BOOL = ISDISABLED()

isImmediate(): BOOL = ISIMMEDIATE() Check if the distribution is equivalent to an Immediate distribution Overrides Distribution.isImmediate(self)

sample(n): X = SAMPLE(N)

class Gamma

Gamma General-purpose continuous distribution with shape and scale

Gamma represents the gamma distribution with shape (alpha) and scale (beta) parameters. This versatile distribution can model various phenomena and includes the exponential and Erlang distributions as special cases. It is commonly used for modeling service times and inter-arrival times.

@brief Gamma distribution with flexible shape and scale parameters

Key characteristics: - Two parameters: shape (alpha) and scale (beta) - Mean = alpha * beta, Variance = alpha * beta² - SCV = 1/alpha (decreases with shape parameter) - Support: (0, ∞) - Includes exponential (alpha=1) and Erlang as special cases

The gamma distribution is used for: - Service time modeling with flexible variability - Waiting time distributions - Inter-arrival time modeling - Reliability and survival analysis - Approximating other continuous distributions

Example: @code service_dist = Gamma(2.0, 1.5); % Shape=2, Scale=1.5, Mean=3.0, SCV=0.5 samples = service_dist.sample(1000); @endcode

Constructor Summary

Gamma(shape, scale)

GAMMA Create a Gamma distribution instance

@brief Creates a Gamma distribution with specified shape and scale @param shape Shape parameter alpha (must be positive) @param scale Scale parameter beta (must be positive) @return self Gamma distribution instance

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T)

evalLST(s): L = EVALLAPLACETRANSFORM(S)

static fitMeanAndSCV(MEAN, SCV): GM = FITMEANANDSCV(MEAN, SCV)

getMean(): EX = GETMEAN()

getProcess(): PROC = GETPROCESS()

getSCV(): SCV = GETSCV()

sample(n): X = SAMPLE(N)

class Exp

Bases: Markovian

Exponential distribution with rate parameter lambda

Fundamental distribution for Markovian queueing systems with memoryless property.

Constructor Summary

Exp(lambda)

EXP Create an exponential distribution instance

@brief Creates an exponential distribution with the given rate parameter @param lambda Rate parameter (must be positive) @return self Exp distribution instance with rate lambda

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T) Evaluate the cumulative distribution function at t AT T

evalLST(s): L = EVALST(S) Evaluate the Laplace-Stieltjes transform of the distribution function at t AT T

static fit(MEAN, SCV, SKEW): EX = FIT(MEAN, SCV, SKEW) Fit the distribution from three standard moments (mean, scv, skewness)

static fitMean(MEAN): EX = FITMEAN(MEAN) Fit exponential distribution with given mean

static fitMeanAndSCV(MEAN, SCV): EX = FITMEANANDSCV(MEAN, SCV) Fit exponential distribution with given mean and squared coefficient of variation (SCV=variance/mean^2)

static fitRate(RATE): EX = FITRATE(RATE) Fit exponential distribution with given rate

static fromMatrix(Lambda): QCELL = FROMMATRIX(LAMBDA) Instantiates a cell array of Exp objects, each with rate given by the entries of the input matrix

getNumberOfPhases(): PHASES = GETNUMBEROFPHASES() Get number of phases in the underpinnning phase-type representation

getSCV()

sample(n): X = SAMPLE(N) Get n samples from the distribution

setMean(MEAN): UPDATEMEAN(SELF,MEAN) Update parameters to match the given mean

setRate(RATE): UPDATERATE(SELF,RATE) Update rate parameter

update(varargin): UPDATE(SELF,VARARGIN) Update parameters to match the first n central moments (n<=4)

class Erlang

Bases: Markovian

Erlang Multi-phase exponential distribution with reduced variability

Erlang represents the distribution of the sum of r independent exponential random variables with the same rate. It models processes with reduced variability compared to exponential distributions and is commonly used for service times with more predictable durations.

@brief Erlang distribution with phase rate and number of phases

Key characteristics: - Sum of r independent exponential phases - Two parameters: phase rate (alpha) and number of phases (r) - Mean = r/alpha, Variance = r/alpha² - SCV = 1/r (always ≤ 1, decreasing with r) - Multi-phase Markovian representation

The Erlang distribution is used for: - Service times with reduced variability - Multi-stage service processes - Modeling deterministic processes approximately - Building blocks for phase-type distributions - Systems where variability decreases with stages

Example: @code service_dist = Erlang(3.0, 4); % 4 phases, rate 3.0 each % Mean = 4/3.0 = 1.33, SCV = 1/4 = 0.25 samples = service_dist.sample(1000); @endcode

Constructor Summary

Erlang(phaseRate, nphases)

ERLANG Create an Erlang distribution instance

@brief Creates an Erlang distribution with specified phase rate and phases @param phaseRate Rate parameter for each exponential phase @param nphases Number of sequential exponential phases (positive integer) @return self Erlang distribution instance

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T)

evalLST(s): L = EVALLAPLACETRANSFORM(S)

static fit(MEAN, SCV, SKEW): ER = FITCENTRAL(MEAN, SCV, SKEW)

static fitMeanAndOrder(MEAN, n): ER = FITMEANANDORDER(MEAN, N)

static fitMeanAndSCV(MEAN, SCV): ER = FITMEANANDSCV(MEAN, SCV)

getMean(): EX = GETMEAN()

getNumberOfPhases(): PHASES = GETNUMBEROFPHASES()

getSCV(): SCV = GETSCV() Get the squared coefficient of variation of the distribution (SCV = variance / mean^2)

setMean(MEAN): UPDATEMEAN(SELF,MEAN) Update parameters to match the given mean

class EmpiricalCDF

Bases: Distribution

Empirical Cdf for a distribution

Constructor Summary

EmpiricalCDF(xdata, cdfdata): SELF = EMPIRICAL(data)

Property Summary

data

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T) Evaluate the cumulative distribution function at t

evalLST(s): L = EVALLST(S) Evaluate the Laplace transform of the distribution function at t

getMean(): MEAN = GETMEAN() Get distribution mean

getMoments()

getRate(): RATE = GETRATE() Get distribution rate

getSCV(): SCV = GETSCV() Get distribution squared coefficient of variation (SCV = variance / mean^2)

getSkewness(): SKEW = GETSKEWNESS() Get distribution skewness

getVar(): VAR = GETVAR() Get distribution variance

sample(n): X = SAMPLE(n) Get n samples from the distribution

class DiscreteDistribution

Bases: Distribution

An abstract class for continuous distributions

class Disabled

Bases: ContinuousDistribution, DiscreteDistribution

A distribution that is not configured

Constructor Summary

Disabled(): SELF = DISABLED()

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T)

evalPMF(k): P = EVALPMF(K)

static getInstance()

getMean(): EX = GETMEAN()

getSCV(): SCV = GETSCV()

isContinuous(): BOOL = ISCONTINUOUS()

isDisabled(): BOOL = ISDISABLED()

isDiscrete(): BOOL = ISDISCRETE()

sample(n): X = SAMPLE(N)

class Det

Bases: ContinuousDistribution, DiscreteDistribution

Det Deterministic distribution with constant value

Det represents a deterministic distribution that always produces the same constant value. This distribution has zero variance and is commonly used for modeling constant service times, fixed delays, or deterministic processing requirements in queueing systems.

@brief Deterministic distribution with constant value and zero variance

Key characteristics: - Single parameter: constant value t - Zero variance (SCV = 0) - Mean = Variance = t - All samples equal to t - Degenerate case of all distributions

The deterministic distribution is used for: - Constant service times - Fixed processing delays - Deterministic timeouts - Modeling perfectly predictable processes - Lower bound analysis in queueing systems

Example: @code constant_service = Det(2.5); % Always takes exactly 2.5 time units samples = constant_service.sample(100); % All samples = 2.5 @endcode

Constructor Summary

Det(t)

DET Create a deterministic distribution instance

@brief Creates a deterministic distribution with constant value t @param t Constant value (must be non-negative for time-based processes) @return self Det distribution instance with constant value t

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T)

evalLST(s): L = EVALST(S)

static fitMean(t)

getMean(): EX = GETMEAN()

getProcess(): PROC = GETPROCESS()

getSCV(): SCV = GETSCV()

sample(n): X = SAMPLE(N)

class ContinuousDistribution

Bases: Distribution

ContinuousDistribution Abstract base class for continuous distributions

ContinuousDistribution provides the common interface and functionality for all continuous-valued statistical distributions. It extends the base Distribution class with methods specific to continuous random variables such as Laplace-Stieltjes transform evaluation.

@brief Abstract base class for continuous-valued distributions

Key characteristics: - Abstract interface for continuous distributions - Support over continuous intervals (often (0,∞) or ℝ) - Provides Laplace-Stieltjes transform interface - Foundation for exponential, gamma, uniform, etc. - Integrates with queueing theory analysis methods

Method Summary

evalLST(s)

class Binomial

A Binomial distribution

Constructor Summary

Binomial(n, p): SELF = BINOMIAL(N, P)

Method Summary

evalCDF(t): FT = EVALCDF(SELF,T)

evalLST(s): L = EVALST(S) Evaluate the Laplace-Stieltjes transform of the distribution function at s For Binomial(n, p), LST(s) = (1 - p + p*e^(-s))^n

evalPMF(k): P = EVALPMF(K)

getMean(): EX = GETMEAN()

getProcess(): PROC = GETPROCESS()

getSCV(): SCV = GETSCV()

getVar(): V = GETVAR()

sample(nsamples): X = SAMPLE(N)

class Bernoulli