Matrix-Analytic Methods
QBD processes and matrix-analytic solutions.
The mam module implements matrix-analytic methods for analyzing queues with
structured Markov chains, including QBD, MAP, RAP, and G/M/1 and M/G/1 type processes.
Key function categories:
MAP analysis:
map_pie(),map_mean(),map_var(),map_scv(),map_skew()MAP fitting:
map2_fit(),mmpp2_fit(),aph_fit(),aph2_fit()PH distributions: Phase-type analysis and fitting
Transformations:
map_scale(),map_normalize(),map_timereverse()QBD methods:
qbd_R(),qbd_mapmap1(),qbd_raprap1()Compression:
compress_adaptive(),compress_spectral()
Matrix-Analytic Methods (MAM) for MAP/PH distributions.
Native Python implementations for analyzing Markovian Arrival Processes (MAPs), Phase-Type (PH) distributions, and related matrix-analytic methods.
- Key algorithms:
map_piq: CTMC steady-state of MAP map_pie: Embedded DTMC steady-state map_lambda: Arrival rate computation map_mean, map_var, map_scv: Moment computations solver_mam_map_bmap_1: MAP/BMAP/1 queue solver using GI/M/1-type ETAQA solver_mam_bmap_map_1: BMAP/MAP/1 queue solver using M/G/1-type ETAQA
- map_infgen(D0, D1)[source]
Compute the infinitesimal generator of a MAP.
The generator Q = D0 + D1 represents the underlying CTMC.
- map_piq(D0, D1=None)[source]
Compute steady-state distribution of the underlying CTMC of a MAP.
Solves πQ = 0 where Q = D0 + D1 is the generator.
- map_pie(D0, D1=None)[source]
Compute equilibrium distribution of embedded DTMC.
The embedded DTMC has transition matrix P = (-D0)^{-1} * D1. Its steady-state is π_e = π * D1 / (π * D1 * e).
- map_lambda(D0, D1=None)[source]
Compute the arrival rate (λ) of a MAP.
The arrival rate is λ = π * D1 * e where π is the steady-state and e is the column vector of ones.
- map_mean(D0, D1=None)[source]
Compute mean inter-arrival time of a MAP.
The mean is 1/λ where λ is the arrival rate.
- map_var(D0, D1=None)[source]
Compute variance of inter-arrival times of a MAP.
Var[X] = E[X²] - E[X]²
Uses map_moment for consistency with Kotlin implementation.
- map_scv(D0, D1=None)[source]
Compute squared coefficient of variation (SCV) of a MAP.
SCV = Var[X] / E[X]² = (E[X²] - E[X]²) / E[X]²
- map_moment(D0, D1, k)[source]
Compute the k-th moment of inter-arrival time distribution.
E[X^k] = k! * π_e * (-D0)^{-k} * e
where π_e is the embedded DTMC steady-state.
- map_scale(D0, D1, factor)[source]
Scale a MAP by a given factor.
Scaling changes the time scale: D0’ = factor * D0, D1’ = factor * D1.
- map_isfeasible(D0, D1, tolerance=1e-10)[source]
Check if (D0, D1) form a valid MAP.
A valid MAP requires: - D0 has non-positive diagonal and non-negative off-diagonal - D1 has non-negative elements - D0 + D1 is a valid generator (row sums = 0)
- map_exponential(mean)[source]
Create a MAP representation of an exponential distribution (MATLAB-style).
This is a MATLAB-compatible wrapper that takes mean instead of rate.
- Parameters:
mean (float) – Mean inter-arrival time (= 1/λ)
- Returns:
Tuple of (D0, D1) matrices representing Exp(1/mean)
- Return type:
Examples
>>> D0, D1 = map_exponential(2) # Poisson process with rate λ=0.5
- map_erlang(mean, k)[source]
Create a MAP representation of an Erlang-k distribution (MATLAB-style).
This is a MATLAB-compatible wrapper that takes mean as first argument.
- Parameters:
- Returns:
Tuple of (D0, D1) matrices representing Erlang-k with given mean
- Return type:
Examples
>>> D0, D1 = map_erlang(2, 3) # Erlang-3 with mean 2
- map_hyperexp(probs, means)[source]
Create a MAP representation of a hyperexponential distribution (MATLAB-style).
This is a MATLAB-compatible wrapper that takes means instead of rates.
- map_gamma(mean, n)[source]
Create a MAP approximation of a Gamma distribution.
Uses Erlang-n with matching mean as an approximation.
- map_sumind(maps)[source]
Compute the sum of independent MAPs.
Creates a MAP representing the sum (concatenation) of independent random variables represented by the input MAPs.
- Parameters:
maps (list) – List of MAPs, each as (D0, D1) tuple
- Returns:
Tuple of (D0, D1) matrices representing the sum
- Return type:
Examples
>>> # Sum of exponential and Erlang-2 >>> MAP1 = map_exponential(1.0) >>> MAP2 = map_erlang(1.0, 2) >>> D0, D1 = map_sumind([MAP1, MAP2])
- map_cdf(D0, D1, points)[source]
Compute cumulative distribution function of inter-arrival times.
F(t) = 1 - π_e * exp(D0*t) * e
- Parameters:
- Returns:
CDF values at specified points
- Return type:
Examples
>>> map_cdf(D0, D1, 1.0) # Returns P(T <= 1) >>> map_cdf(D0, D1, [1.0, 5.0]) # Returns [P(T<=1), P(T<=5)]
- map_pdf(D0, D1, points)[source]
Compute probability density function of inter-arrival times.
f(t) = π_e * exp(D0*t) * (-D0) * e
- Parameters:
- Returns:
PDF values at specified points
- Return type:
Examples
>>> map_pdf(D0, D1, [0.5, 1.0, 2.0])
- map_sample(D0, D1, n_samples, rng=None)[source]
Generate random samples from a MAP distribution.
- Parameters:
D0 (ndarray) – Hidden transition matrix
D1 (ndarray) – Visible transition matrix
n_samples (int) – Number of samples to generate
rng (numpy.random.Generator | None) – Optional random number generator
- Returns:
Array of inter-arrival times
- Return type:
- map_acf(D0, D1, lags=1)[source]
Compute autocorrelation coefficients of inter-arrival times.
- Parameters:
- Returns:
Array of autocorrelation coefficients at specified lags
- Return type:
Examples
>>> map_acf(D0, D1) # lag-1 autocorrelation >>> map_acf(D0, D1, np.arange(1, 11)) # first 10 autocorrelations
- map_idc(D0, D1)[source]
Compute the asymptotic index of dispersion.
I = SCV * (1 + 2 * sum_{k=1}^{inf} rho_k)
where SCV is the squared coefficient of variation and rho_k is the lag-k autocorrelation coefficient.
- map_count_var(D0, D1, t)[source]
Compute variance of counting process at resolution t.
- Parameters:
- Returns:
Variance of arrivals in (0, t]
- Return type:
- Reference:
He and Neuts, “Markov chains with marked transitions”, 1998
- map_varcount(D0, D1, tset)[source]
Compute variance of counting process (alternative implementation).
- map_count_moment(D0, D1, t, orders)[source]
Compute power moments of counts at resolution t.
Uses numerical differentiation of the moment generating function.
- map_mmpp2(mean, scv, skew=-1, acf1=-1)[source]
Fit an MMPP(2) as a MAP.
- Parameters:
- Returns:
Tuple of (D0, D1) matrices
- Return type:
Examples
>>> D0, D1 = map_mmpp2(1, 2, -1, 0.2) # Minimal skewness, ACF=0.2 >>> D0, D1 = map_mmpp2(1, 2, -1, -1) # Minimal skewness, max ACF
- map_gamma2(D0, D1)[source]
Compute the second largest eigenvalue of embedded DTMC.
This is the autocorrelation decay rate.
- map_randn(k, mu=(1.0, 1.0), sigma=(0.5, 0.5))[source]
Generate a random MAP with normally distributed elements.
- map_pntiter(D0, D1, na, t, M=None)[source]
Compute probability of na arrivals in interval [0, t].
Uses iterative bisection method (Neuts and Li).
- map2_fit(e1, e2, e3=-1.0, g2=0.0)[source]
Fit a MAP(2) distribution to moments and autocorrelation.
Based on: A. Heindl, G. Horvath, K. Gross “Explicit inverse characterization of acyclic MAPs of second order”
- Parameters:
- Returns:
MAP: Tuple (D0, D1) if successful, None if failed
error_code: 0 for success, >0 for various errors
- Return type:
Tuple of (MAP, error_code) where
- Error codes:
0: Success 10: Mean out of bounds 20: Correlated exponential 30: h2 out of bounds 40: h3 out of bounds 51-54: g2 out of bounds
- map_joint(D0, D1, a, i)[source]
Compute joint moments of a MAP.
E[(X_{a1})^{i1} * (X_{a1+a2})^{i2} * …]
- map_issym(D0, D1=None)[source]
Check if MAP contains symbolic elements.
In Python, this always returns False as we use numpy arrays.
- map_feastol()[source]
Get tolerance for feasibility checks.
- Returns:
Tolerance exponent (10^(-feastol))
- Return type:
- map_largemap()[source]
Get threshold for “large” MAP where exact computation is expensive.
- Returns:
- Return type:
Order threshold (default
- mmap_normalize(D0, D_list)[source]
Normalize an MMAP to ensure -diag(D0) = rowsums(D1 + D2 + … + Dm).
This ensures D0 represents true hidden transitions (negative diagonal).
- mmap_super_safe(mmap_list, maxorder=1000, method='default')[source]
Safely superpose multiple MMAPs into a single MMAP.
Combines multiple MMAPs while ensuring the order of the resulting MMAP does not exceed maxorder. If combining would exceed maxorder, alternative methods are used (fitting to exponential or simplified MAP).
- Parameters:
- Returns:
Tuple of (D0_super, D_list_super) for the superposed MMAP
- Return type:
- Algorithm:
Sort MMAPs by squared coefficient of variation (SCV) of unmarked process
Iteratively combine: start with simplest, add others in order
If Kronecker product would exceed maxorder, use simpler approximations
- mmap_mark(D0, D_list, prob)[source]
Reclassify arrivals in an MMAP according to a probability matrix.
Converts an MMAP with K classes to one with R classes based on a KxR probability matrix that describes how arrivals are reclassified.
- Parameters:
- Returns:
Tuple of (D0_new, D_list_new) where D_list_new has R matrices
- Return type:
Example
If prob = [[0.8, 0.2], [0.3, 0.7]], a 2-class MMAP becomes a 2-class MMAP where 80% of class-1 arrivals are now class-1, 20% become class-2, etc.
- mmap_scale(D0, D_list, M, max_iter=30)[source]
Scale the mean inter-arrival times of an MMAP.
Adjusts all matrices to achieve specified mean inter-arrival times.
- Parameters:
- Returns:
Tuple of (D0_scaled, D_list_scaled)
- Return type:
- Algorithm:
For single M: Scale all matrices uniformly by ratio = M_old / M For vector M: Use iterative coordinate descent to find scaling factors
- mmap_hide(D0, D_list, types)[source]
Hide (remove) specified arrival classes from an MMAP.
The hidden classes are set to zero matrices, effectively removing them from observation while maintaining the underlying stochastic process.
- Parameters:
- Returns:
Tuple of (D0_hidden, D_list_hidden) where specified classes are zero
- Return type:
Example
mmap_hide(D0, [D1, D2, D3], types=[1]) hides class 2 (index 1)
- mmap_compress(D0, D_list, method='default', target_order=None)[source]
Compress an MMAP using various approximation methods.
Reduces the state space of an MMAP while preserving key statistical properties (moments, inter-arrival time distribution, etc.).
- Parameters:
- Returns:
Tuple of (D0_compressed, D_list_compressed)
- Return type:
Notes
The “default” mixture method fits each class independently to MMPP(2) if possible, then combines them. This is fast and often effective.
More advanced methods (m3pp, mamap2) can be added as extensions.
- mmap_exponential(lambda_rate, nclasses=None)[source]
Create an exponential MMAP (Poisson arrival process).
- mmap_issym(mmap)[source]
Check if an MMAP contains symbolic elements.
In Python/NumPy, we don’t have built-in symbolic support like MATLAB, so this always returns False for numpy arrays.
- mmap_isfeasible(mmap, tol=None)[source]
Check whether an MMAP is feasible up to the given tolerance.
Checks: - Elements are real (no imaginary parts) - D0 + D1 rows sum to zero (generator property) - Diagonal of D0 is negative - Off-diagonal of D0 is non-negative - All D1c matrices are non-negative - D1 = D11 + D12 + … + D1C
- mmap_count_lambda(mmap)[source]
Compute the arrival rate of the counting process for each class.
The rate for class k is λ_k = π * D_k * e where π is the steady-state of the underlying CTMC.
- mmap_pie(mmap)[source]
Compute the stationary probability of the DTMC embedded at restart instants after an arrival of each class.
- mmap_embedded(mmap)[source]
Compute the embedded DTMC transition matrices for each class.
For class k, the embedded matrix is Pc_k = (-D0)^{-1} * D_{2+k}
- mmap_sample(mmap, n_samples, pi=None, seed=None)[source]
Generate samples from a Marked MAP.
- Parameters:
- Returns:
T: Array of inter-arrival times
A: Array of class labels (1-indexed for MATLAB compatibility)
- Return type:
Tuple of (T, A) where
- mmap_sum(mmap, n)[source]
Create an MMAP representing the sum of n identical MMAPs (n-fold sum).
This creates an MMAP whose inter-arrival times are distributed as the sum of n inter-arrival times from the original MMAP.
- mmap_super(mmap_a, mmap_b=None, opt='default')[source]
Superpose two or more MMAPs.
- Parameters:
- Returns:
Superposed MMAP
- Return type:
- mmap_mixture(alpha, maps)[source]
Create a probabilistic mixture of MAPs.
Each MAP in the list is selected with probability alpha[i].
- mmap_max(mmap_a, mmap_b, k)[source]
Create MMAP for the maximum of arrivals from two synchronized MMAPs.
This models a synchronization queue where arrivals from both sources must be paired before release.
- mmap_maps(mmap)[source]
Extract K MAPs, one for each class of the MMAP[K] process.
For class k, the MAP is {D0 + D1 - D_{2+k}, D_{2+k}}
- mmap_count_idc(mmap, t)[source]
Compute the per-class Index of Dispersion of Counts at resolution t.
IDC = Var[N(t)] / E[N(t)]
- mmap_count_mcov(mmap, t)[source]
Compute the count covariance between each pair of classes at time scale t.
- mmap_idc(mmap)[source]
Compute the asymptotic Index of Dispersion of Counts for each class.
This is the limit of IDC(t) as t -> infinity.
- mmap_sigma(mmap)[source]
Compute one-step class transition probabilities.
p_{i,j} = P(C_k = j | C_{k-1} = i)
- mmap_sigma2(mmap)[source]
Compute two-step class transition probabilities.
p_{i,j,h} = P(C_k = h | C_{k-1} = j, C_{k-2} = i)
- mmap_forward_moment(mmap, orders, norm=True)[source]
Compute the theoretical forward moments of an MMAP.
Forward moments are E[X^k | previous arrival was class c].
- Parameters:
- Returns:
Matrix of shape (C, len(orders)) where element (c, k) is the order-k forward moment for class c
- Return type:
- mmap_backward_moment(mmap, orders, norm=True)[source]
Compute the theoretical backward moments of an MMAP.
Backward moments are E[X^k | next arrival will be class c].
- Parameters:
- Returns:
Matrix of shape (C, len(orders)) where element (c, k) is the order-k backward moment for class c
- Return type:
- mmap_cross_moment(mmap, k)[source]
Compute the k-th order moment of inter-arrival times between class pairs.
Computes E[X^k | previous = class i, next = class j] for all i, j.
- mmap_modulate(P, HT, MMAP)[source]
Modulate MMAPs in continuous time according to holding time distributions.
Creates an MMAP that switches between different MMAPs according to a continuous-time Markov modulated process.
- ldqbd(Q0, Q1, Q2, options=None)[source]
Solve a level-dependent QBD process.
- Parameters:
Q0 (List[ndarray]) – List of upward transition matrices [Q0^(0), Q0^(1), …, Q0^(N-1)] Q0^(n) has shape (states_n, states_{n+1})
Q1 (List[ndarray]) – List of local transition matrices [Q1^(0), Q1^(1), …, Q1^(N)] Q1^(n) has shape (states_n, states_n)
Q2 (List[ndarray]) – List of downward transition matrices [Q2^(1), Q2^(2), …, Q2^(N)] Q2^(n) has shape (states_n, states_{n-1})
options (LdqbdOptions | None) – LdqbdOptions instance (uses defaults if None)
- Returns:
LdqbdResult containing rate matrices R and stationary distribution pi
- Return type:
- Algorithm:
Compute rate matrices backward using continued fraction recursion
Compute stationary distribution forward using R matrices
- class LdqbdResult(R, pi)[source]
Bases:
objectResult of LDQBD solver.
- R
List of rate matrices R^(1), R^(2), …, R^(N)
- Type:
List[numpy.ndarray]
- pi
Stationary distribution vector [pi_0, pi_1, …, pi_N]
- Type:
- class LdqbdOptions(epsilon=1e-10, max_iter=1000, verbose=False)[source]
Bases:
objectOptions for LDQBD solver.
- class MAPBMAP1Result(mean_queue_length, utilization, mean_response_time, throughput, pi, R, mean_batch_size)[source]
Bases:
objectResult of MAP/BMAP/1 queue analysis.
- solver_mam_map_bmap_1(C0, C1, D)[source]
Solve a MAP/BMAP/1 queue using GI/M/1-type matrix-analytic methods.
The MAP is specified by matrices (C0, C1) where: - C0: transitions without arrivals - C1: transitions triggering arrivals
The BMAP for service is specified by matrices {D0, D1, D2, …, DK} where: - D0: transitions without service completions - Dk: transitions triggering batch service of k customers (k >= 1)
The GI/M/1-type structure for MAP/BMAP/1 is: ```
B1 A0 0 0 … B2 A1 A0 0 …
- Q = B3 A2 A1 A0 …
…
- Where:
A0 = C1 otimes I_ms (MAP arrival, level +1) A1 = C0 otimes I_ms + I_ma otimes D0 (phase changes, level 0) A_{k+1} = I_ma otimes D_k (batch size k service, level -k)
- Parameters:
- Returns:
MAPBMAP1Result with performance metrics
- Return type:
- class BMAPMAP1Result(mean_queue_length, utilization, mean_response_time, throughput, pi, G, mean_batch_size)[source]
Bases:
objectResult of BMAP/MAP/1 queue analysis.
- solver_mam_bmap_map_1(D, S0, S1)[source]
Solve a BMAP/MAP/1 queue using M/G/1-type matrix-analytic methods.
The BMAP is specified by matrices {D0, D1, D2, …, DK} where: - D0: transitions without arrivals - Dk: transitions triggering batch size k arrivals (k >= 1)
The MAP for service is specified by matrices (S0, S1) where: - S0: transitions without service completions - S1: transitions triggering service completions
The M/G/1-type structure for BMAP/MAP/1 is: ```
A0 = I_ma otimes S1 (service completion, level -1) A1 = D0 otimes I_ms + I_ma otimes S0 (phase changes, level 0) A_{k+1} = D_k otimes I_ms (batch arrival size k, level +k)
- Parameters:
- Returns:
BMAPMAP1Result with performance metrics
- Return type:
- qbd_R(B, L, F, iter_max=100000, tol=1e-12)[source]
Compute QBD rate matrix R using successive substitutions.
- Solves the matrix quadratic equation:
R^2 * A_{-1} + R * A_0 + A_1 = 0
where A_{-1} = B, A_0 = L, A_1 = F.
- Parameters:
- Returns:
Rate matrix R
- Return type:
References
Original MATLAB: matlab/src/api/mam/qbd_R.m
- qbd_R_logred(B, L, F, iter_max=1000, tol=1e-14)[source]
Compute QBD rate matrix R using logarithmic reduction.
Uses the logarithmic reduction algorithm which has quadratic convergence compared to linear convergence of successive substitutions.
- Parameters:
- Returns:
Rate matrix R
- Return type:
References
Original MATLAB: matlab/src/api/mam/qbd_R_logred.m Latouche & Ramaswami, Ch. 8
- qbd_rg(B, L, F, method='logred', iter_max=1000, tol=1e-14)[source]
Compute both R and G matrices for a QBD process.
- G is the minimal non-negative solution to:
A_1 * G^2 + A_0 * G + A_{-1} = 0
- R is the minimal non-negative solution to:
R^2 * A_{-1} + R * A_0 + A_1 = 0
- Parameters:
- Returns:
QBDResult with R, G, U, and eta (caudal characteristic)
- Return type:
References
Original MATLAB: matlab/src/api/mam/qbd_rg.m
- qbd_blocks_mapmap1(D0_arr, D1_arr, D0_srv, D1_srv)[source]
Construct QBD blocks for a MAP/MAP/1 queue.
Builds the backward (B), local (L), and forward (F) transition blocks for the QBD representation of a MAP/MAP/1 queue.
- Parameters:
- Returns:
Tuple of (B, L, F) QBD blocks
- Return type:
References
Original MATLAB: matlab/src/api/mam/qbd_mapmap1.m
- qbd_bmapbmap1(MAPa, pbatcha, MAPs)[source]
Compute QBD blocks for a BMAP/BMAP/1 queue.
Constructs the QBD (Quasi-Birth-Death) transition blocks for a BMAP/BMAP/1 queue with batch arrivals.
- Parameters:
- Returns:
A0: Local transition block A_1: Downward transition block A1_list: List of upward transition blocks for each batch size B0: Initial boundary local block B1_list: List of boundary upward blocks for each batch size
- Return type:
Tuple of (A0, A_1, A1_list, B0, B1_list) where
References
Original MATLAB: matlab/src/api/mam/qbd_bmapbmap1.m
- qbd_mapmap1(MAPa, MAPs, util=None)[source]
Analyze a MAP/MAP/1 queue using QBD methods.
Solves a MAP/MAP/1 queue using Quasi-Birth-Death process methods, computing throughput, queue length, utilization, and other metrics.
- Parameters:
- Returns:
Tuple of (XN, QN, UN, pqueue, R, eta, G, A_1, A0, A1, U, MAPs_scaled) where:
XN: System throughput QN: Mean queue length UN: Utilization pqueue: Queue length distribution R: Rate matrix R eta: Caudal characteristic (spectral radius of R) G: Rate matrix G A_1: Downward transition block A0: Local transition block A1: Upward transition block U: Matrix U MAPs_scaled: Scaled service process
- Return type:
Tuple[float, float, float, ndarray, ndarray, float | None, ndarray, ndarray, ndarray, ndarray, ndarray, Tuple[ndarray, ndarray]]
References
Original MATLAB: matlab/src/api/mam/qbd_mapmap1.m
- qbd_raprap1(RAPa, RAPs, util=None)[source]
Analyze a RAP/RAP/1 queue using QBD methods.
Solves a RAP/RAP/1 queue (Rational Arrival Process) using QBD methods, computing throughput, queue length, utilization, and other metrics.
Note: This is a simplified implementation. For full RAP analysis, specialized algorithms like those in SMCSolver are needed.
- Parameters:
- Returns:
XN: System throughput QN: Mean queue length UN: Utilization pqueue: Queue length distribution R: Rate matrix R eta: Caudal characteristic G: Rate matrix G B: Backward transition block L: Local transition block F: Forward transition block
- Return type:
Tuple of (XN, QN, UN, pqueue, R, eta, G, B, L, F) where
References
Original MATLAB: matlab/src/api/mam/qbd_raprap1.m
- qbd_setupdelayoff(lambda_val, mu, alpharate, alphascv, betarate, betascv)[source]
Analyze queue with setup delay and turn-off phases.
Performs queue-length analysis for a queueing system with setup delay (warm-up) and turn-off periods using QBD methods.
The system operates as follows: 1. When empty and job arrives, server enters setup phase 2. After setup, server becomes active and serves jobs 3. When queue empties, server enters turn-off phase 4. After turn-off, server becomes idle
- Parameters:
- Returns:
Average queue length QN
- Return type:
References
Original MATLAB: matlab/src/api/mam/qbd_setupdelayoff.m
- map_ccdf_derivative(MAP, i)[source]
Compute derivative at zero of a MAP’s complementary CDF.
Calculates the i-th derivative at t=0 of the complementary cumulative distribution function (CCDF) of the inter-arrival time distribution.
Formula: ν_i = π_e * D_0^i * e
where π_e is the embedded stationary vector and e is the column vector of ones.
- Parameters:
- Returns:
Value of the i-th derivative at zero
- Return type:
References
Original MATLAB: matlab/src/api/mam/map_ccdf_derivative.m A. Horvath et al., “A Joint Moments Based Analysis of Networks of MAP/MAP/1 Queues”
- map_jointpdf_derivative(MAP, iset)[source]
Compute partial derivative at zero of a MAP’s joint PDF.
Calculates the partial derivative at t=0 of the joint probability density function of consecutive inter-arrival times.
For index set {i_1, i_2, …, i_k}: γ = π_e * D_0^{i_1} * D_1 * D_0^{i_2} * D_1 * … * D_0^{i_k} * D_1 * e
- Parameters:
- Returns:
Value of the partial derivative at zero
- Return type:
References
Original MATLAB: matlab/src/api/mam/map_jointpdf_derivative.m A. Horvath et al., “A Joint Moments Based Analysis of Networks of MAP/MAP/1 Queues”
- map_factorial_moment(MAP, k)[source]
Compute the k-th factorial moment of a MAP.
The k-th factorial moment is computed using derivatives of the inter-arrival time distribution.
- Parameters:
- Returns:
k-th factorial moment
- Return type:
References
Based on MAP moment formulas from matrix-analytic methods
- map_joint_moment(MAP, k, l)[source]
Compute the (k,l)-th joint moment of consecutive inter-arrival times.
E[X_n^k * X_{n+1}^l] for a MAP with inter-arrival times X_n.
- Parameters:
- Returns:
Joint moment E[X_n^k * X_{n+1}^l]
- Return type:
References
Based on joint moment formulas for MAPs
- map_m1ps_cdf_respt(C, D, mu, x, epsilon=1e-11, epsilon_prime=1e-10)[source]
Compute complementary sojourn time CDF for MAP/M/1-PS queue.
- Based on:
Masuyama, H., & Takine, T. (2003). Sojourn time distribution in a MAP/M/1 processor-sharing queue.
- Parameters:
- Returns:
Complementary CDF values W_bar(x) = Pr[W > x] at each point in x.
- Return type:
Example
>>> # Exponential arrivals (Poisson process with rate lambda=0.5) >>> C = np.array([[-0.5]]) >>> D = np.array([[0.5]]) >>> mu = 1.0 >>> x = np.array([0.0, 0.5, 1.0, 2.0, 5.0]) >>> cdf = map_m1ps_cdf_respt(C, D, mu, x)
- map_compute_R(C, D, mu)[source]
Compute rate matrix R for MAP/M/1 queue.
- Computes the minimal nonnegative solution of the matrix equation:
D + R(C - mu*I) + mu*R^2 = 0
This matrix is used in the analysis of MAP/M/1 queues based on quasi-birth-death processes.
- Parameters:
- Returns:
M x M rate matrix (minimal nonnegative solution)
- Return type:
References
Original MATLAB: matlab/src/api/map/map_compute_R.m Masuyama, H., & Takine, T. (2003). Sojourn time distribution in a MAP/M/1 processor-sharing queue.
- map_m1ps_h_recursive(C, D, mu, N, K)[source]
Recursive computation of h_{n,k} coefficients for MAP/M/1-PS.
- The h_{n,k} vectors satisfy the recursion (Theorem 1 in paper):
h_{n,0} = e (vector of ones), for n = 0, 1, … h_{n,k+1} = 1/(theta+mu) * [n*mu/(n+1) * h_{n-1,k} + (theta*I + C) * h_{n,k}
D * h_{n+1,k}]
where h_{-1,k} = 0 for all k
- Parameters:
C (ndarray | list) – M x M matrix governing MAP transitions without arrivals
D (ndarray | list) – M x M matrix governing MAP transitions with arrivals
mu (float) – Service rate (scalar)
N (int) – Maximum value of n to compute (determines rows)
K (int) – Maximum value of k to compute (determines columns)
- Returns:
2D list h[n][k] of (M,1) arrays containing h_{n,k}
- Return type:
References
Original MATLAB: matlab/src/api/map/map_m1ps_h_recursive.m Masuyama, H., & Takine, T. (2003).
- map_m1ps_sojourn(C, D, mu, x, epsilon=1e-11, epsilon_prime=1e-10)[source]
Compute sojourn time distribution in MAP/M/1-PS queue.
Alias for map_m1ps_cdf_respt.
- Parameters:
- Returns:
Complementary CDF values W_bar(x) = Pr[W > x]
- Return type:
References
Original MATLAB: matlab/src/api/map/map_m1ps_sojourn.m
- mmdp_isfeasible(Q, R, tol=1e-10)[source]
Check if (Q, R) defines a valid MMDP.
Requirements: - Q must be a valid generator (square, row sums = 0, proper signs) - R must be diagonal with non-negative entries - Q and R must have compatible dimensions
- Parameters:
- Returns:
True if (Q, R) defines a valid MMDP, False otherwise
- Return type:
Examples
>>> Q = np.array([[-0.5, 0.5], [0.3, -0.3]]) >>> R = np.array([[2.0, 0], [0, 5.0]]) >>> mmdp_isfeasible(Q, R) True