"""
Discrete probability distributions for LINE (pure Python).
This module provides discrete distribution implementations including
Poisson, geometric, binomial, and other common distributions.
"""
from typing import Optional, Tuple, Union
import numpy as np
from scipy.special import factorial
class _LazyStats:
def __getattr__(self, name):
from scipy import stats as _scipy_stats
globals()['stats'] = _scipy_stats
return getattr(_scipy_stats, name)
stats = _LazyStats()
from .base import DiscreteDistribution
from line_solver.api.io.errors import LineInputError
[docs]
class Poisson(DiscreteDistribution):
"""
Poisson distribution.
The Poisson distribution models the number of events occurring
in a fixed interval of time or space.
Args:
lambda_: Rate parameter (mean = variance = lambda).
"""
def __init__(self, lambda_: float):
super().__init__()
self._name = 'Poisson'
if lambda_ <= 0:
raise ValueError("Lambda must be positive")
self._lambda = lambda_
@property
def lambda_(self) -> float:
"""Get the rate parameter."""
return self._lambda
[docs]
def getMean(self) -> float:
"""Get the mean (equals lambda)."""
return self._lambda
[docs]
def getVar(self) -> float:
"""Get the variance (equals lambda)."""
return self._lambda
[docs]
def getSCV(self) -> float:
"""Get the SCV (1/lambda)."""
return 1.0 / self._lambda
[docs]
def getSkew(self) -> float:
"""Get the skewness."""
return 1.0 / np.sqrt(self._lambda)
[docs]
def getSupport(self) -> Tuple[float, float]:
"""Get the support [0, inf)."""
return (0.0, float('inf'))
[docs]
def evalPMF(self, x: int) -> float:
"""Evaluate the probability mass function at point x."""
if x < 0:
return 0.0
return stats.poisson.pmf(x, self._lambda)
[docs]
def evalCDF(self, x: float) -> float:
"""Evaluate the CDF at point x."""
if x < 0:
return 0.0
return stats.poisson.cdf(int(x), self._lambda)
[docs]
def sample(self, n: int = 1, rng: Optional[np.random.Generator] = None) -> np.ndarray:
"""Generate random samples."""
if rng is None:
rng = np.random.default_rng()
return rng.poisson(lam=self._lambda, size=n)
[docs]
class Geometric(DiscreteDistribution):
"""
Geometric distribution.
The geometric distribution models the number of trials up to and
including the first success in a sequence of Bernoulli trials, so its
support is {1, 2, ...}.
Args:
p: Probability of success on each trial (0 < p <= 1).
"""
def __init__(self, p: float):
super().__init__()
self._name = 'Geometric'
if p <= 0 or p > 1:
raise ValueError("Probability p must be in (0, 1]")
self._p = p
@property
def p(self) -> float:
"""Get the probability of success."""
return self._p
[docs]
def getMean(self) -> float:
"""Get the mean (1/p for number of trials)."""
return 1.0 / self._p
[docs]
def getVar(self) -> float:
"""Get the variance."""
return (1 - self._p) / (self._p ** 2)
[docs]
def getSkew(self) -> float:
"""Get the skewness."""
return (2 - self._p) / np.sqrt(1 - self._p)
[docs]
def getSupport(self) -> Tuple[float, float]:
"""Get the support [1, inf)."""
return (1.0, float('inf'))
[docs]
def evalPMF(self, x: int) -> float:
"""Evaluate the PMF at point x (number of trials until first success)."""
if x < 1:
return 0.0
return self._p * (1 - self._p) ** (x - 1)
[docs]
def evalCDF(self, x: float) -> float:
"""Evaluate the CDF at point x."""
if x < 1:
return 0.0
return 1 - (1 - self._p) ** int(x)
[docs]
def sample(self, n: int = 1, rng: Optional[np.random.Generator] = None) -> np.ndarray:
"""Generate random samples."""
if rng is None:
rng = np.random.default_rng()
return rng.geometric(p=self._p, size=n)
[docs]
class Binomial(DiscreteDistribution):
"""
Binomial distribution.
The binomial distribution models the number of successes in
n independent Bernoulli trials.
Args:
n: Number of trials.
p: Probability of success on each trial.
"""
def __init__(self, n: int, p: float):
super().__init__()
self._name = 'Binomial'
if n < 0:
raise ValueError("Number of trials must be non-negative")
if p < 0 or p > 1:
raise ValueError("Probability p must be in [0, 1]")
self._n = n
self._p = p
@property
def n(self) -> int:
"""Get the number of trials."""
return self._n
@property
def p(self) -> float:
"""Get the probability of success."""
return self._p
[docs]
def getMean(self) -> float:
"""Get the mean (n * p)."""
return self._n * self._p
[docs]
def getVar(self) -> float:
"""Get the variance (n * p * (1 - p))."""
return self._n * self._p * (1 - self._p)
[docs]
def getSkew(self) -> float:
"""Get the skewness."""
if self._n == 0 or self._p == 0 or self._p == 1:
return 0.0
return (1 - 2 * self._p) / np.sqrt(self._n * self._p * (1 - self._p))
[docs]
def getSupport(self) -> Tuple[float, float]:
"""Get the support [0, n]."""
return (0.0, float(self._n))
[docs]
def evalPMF(self, x: int) -> float:
"""Evaluate the PMF at point x."""
if x < 0 or x > self._n:
return 0.0
return stats.binom.pmf(x, self._n, self._p)
[docs]
def evalCDF(self, x: float) -> float:
"""Evaluate the CDF at point x."""
if x < 0:
return 0.0
if x >= self._n:
return 1.0
return stats.binom.cdf(int(x), self._n, self._p)
[docs]
def sample(self, n: int = 1, rng: Optional[np.random.Generator] = None) -> np.ndarray:
"""Generate random samples."""
if rng is None:
rng = np.random.default_rng()
return rng.binomial(n=self._n, p=self._p, size=n)
class NegBinomial(DiscreteDistribution):
"""
Negative Binomial distribution.
Models the number of failures before r successes in
a sequence of Bernoulli trials.
Args:
r: Number of successes required.
p: Probability of success on each trial.
"""
def __init__(self, r: int, p: float):
super().__init__()
self._name = 'NegBinomial'
if r < 1:
raise ValueError("Number of successes r must be at least 1")
if p <= 0 or p > 1:
raise ValueError("Probability p must be in (0, 1]")
self._r = r
self._p = p
@property
def r(self) -> int:
"""Get the number of successes required."""
return self._r
@property
def p(self) -> float:
"""Get the probability of success."""
return self._p
def getMean(self) -> float:
"""Get the mean (r * (1 - p) / p)."""
return self._r * (1 - self._p) / self._p
def getVar(self) -> float:
"""Get the variance."""
return self._r * (1 - self._p) / (self._p ** 2)
def getSkew(self) -> float:
"""Get the skewness."""
return (2 - self._p) / np.sqrt(self._r * (1 - self._p))
def getSupport(self) -> Tuple[float, float]:
"""Get the support [0, inf)."""
return (0.0, float('inf'))
def evalPMF(self, x: int) -> float:
"""Evaluate the PMF at point x."""
if x < 0:
return 0.0
return stats.nbinom.pmf(x, self._r, self._p)
def evalCDF(self, x: float) -> float:
"""Evaluate the CDF at point x."""
if x < 0:
return 0.0
return stats.nbinom.cdf(int(x), self._r, self._p)
def sample(self, n: int = 1, rng: Optional[np.random.Generator] = None) -> np.ndarray:
"""Generate random samples."""
if rng is None:
rng = np.random.default_rng()
return rng.negative_binomial(n=self._r, p=self._p, size=n)
[docs]
class Zipf(DiscreteDistribution):
"""
Zipf distribution.
The Zipf distribution is a power-law distribution often used
to model rank-frequency relationships (e.g., word frequencies).
Args:
s: Shape parameter (s > 0).
n: Upper bound on support (optional, defaults to large value).
"""
def __init__(self, s: float, n: int = 10000):
super().__init__()
self._name = 'Zipf'
if s <= 0:
raise ValueError("Shape parameter s must be > 0")
if n < 1:
raise ValueError("Upper bound n must be at least 1")
self._s = s
self._n = n
# Compute normalization constant
self._H = np.sum(1.0 / np.arange(1, n + 1) ** s)
@property
def s(self) -> float:
"""Get the shape parameter."""
return self._s
@property
def n(self) -> int:
"""Get the upper bound."""
return self._n
[docs]
def getMean(self) -> float:
"""Get the mean."""
k = np.arange(1, self._n + 1)
return np.sum(k * (1.0 / k ** self._s) / self._H)
[docs]
def getVar(self) -> float:
"""Get the variance."""
k = np.arange(1, self._n + 1)
probs = (1.0 / k ** self._s) / self._H
mean = np.sum(k * probs)
return np.sum(k ** 2 * probs) - mean ** 2
[docs]
def getSupport(self) -> Tuple[float, float]:
"""Get the support [1, n]."""
return (1.0, float(self._n))
[docs]
def evalPMF(self, x: int) -> float:
"""Evaluate the PMF at point x."""
if x < 1 or x > self._n:
return 0.0
return (1.0 / x ** self._s) / self._H
[docs]
def evalCDF(self, x: float) -> float:
"""Evaluate the CDF at point x."""
if x < 1:
return 0.0
if x >= self._n:
return 1.0
k = int(x)
return np.sum(1.0 / np.arange(1, k + 1) ** self._s) / self._H
[docs]
def sample(self, n: int = 1, rng: Optional[np.random.Generator] = None) -> np.ndarray:
"""Generate random samples using inverse transform."""
if rng is None:
rng = np.random.default_rng()
# Build CDF
k = np.arange(1, self._n + 1)
pmf = (1.0 / k ** self._s) / self._H
cdf = np.cumsum(pmf)
# Inverse transform sampling
u = rng.random(size=n)
return np.searchsorted(cdf, u) + 1
class Empirical(DiscreteDistribution):
"""
Empirical discrete distribution from data.
Creates a distribution from observed data by computing
the empirical probability mass function.
Args:
data: Array of observed values.
"""
def __init__(self, data: Union[list, np.ndarray]):
super().__init__()
self._name = 'Empirical'
self._data = np.array(data)
# Compute unique values and their probabilities
unique, counts = np.unique(self._data, return_counts=True)
self._values = unique
self._probs = counts / len(self._data)
self._cdf = np.cumsum(self._probs)
@property
def values(self) -> np.ndarray:
"""Get the unique values."""
return self._values
@property
def probs(self) -> np.ndarray:
"""Get the probabilities for each value."""
return self._probs
def getMean(self) -> float:
"""Get the mean."""
return float(np.sum(self._values * self._probs))
def getVar(self) -> float:
"""Get the variance."""
mean = self.getMean()
return float(np.sum(self._probs * (self._values - mean) ** 2))
def getSupport(self) -> Tuple[float, float]:
"""Get the support [min, max]."""
return (float(np.min(self._values)), float(np.max(self._values)))
def evalPMF(self, x: int) -> float:
"""Evaluate the PMF at point x."""
idx = np.where(self._values == x)[0]
if len(idx) == 0:
return 0.0
return self._probs[idx[0]]
def evalCDF(self, x: float) -> float:
"""Evaluate the CDF at point x."""
idx = np.searchsorted(self._values, x, side='right') - 1
if idx < 0:
return 0.0
return self._cdf[min(idx, len(self._cdf) - 1)]
def sample(self, n: int = 1, rng: Optional[np.random.Generator] = None) -> np.ndarray:
"""Generate random samples."""
if rng is None:
rng = np.random.default_rng()
indices = rng.choice(len(self._values), size=n, p=self._probs)
return self._values[indices]
[docs]
class Bernoulli(DiscreteDistribution):
"""
Bernoulli distribution.
The simplest discrete distribution with two outcomes: success (1)
with probability p and failure (0) with probability 1-p.
Args:
p: Probability of success (0 <= p <= 1).
"""
def __init__(self, p: float):
super().__init__()
self._name = 'Bernoulli'
if p < 0 or p > 1:
raise ValueError("Probability p must be in [0, 1]")
self._p = p
@property
def p(self) -> float:
"""Get the probability of success."""
return self._p
[docs]
def getMean(self) -> float:
"""Get the mean (equals p)."""
return self._p
[docs]
def getVar(self) -> float:
"""Get the variance (p * (1 - p))."""
return self._p * (1 - self._p)
[docs]
def getSkew(self) -> float:
"""Get the skewness."""
if self._p == 0 or self._p == 1:
return 0.0
return (1 - 2 * self._p) / np.sqrt(self._p * (1 - self._p))
[docs]
def getSupport(self) -> Tuple[float, float]:
"""Get the support {0, 1}."""
return (0.0, 1.0)
[docs]
def evalPMF(self, x: int) -> float:
"""Evaluate the PMF at point x."""
if x == 0:
return 1 - self._p
elif x == 1:
return self._p
else:
return 0.0
[docs]
def evalCDF(self, x: float) -> float:
"""Evaluate the CDF at point x."""
if x < 0:
return 0.0
elif x < 1:
return 1 - self._p
else:
return 1.0
[docs]
def sample(self, n: int = 1, rng: Optional[np.random.Generator] = None) -> np.ndarray:
"""Generate random samples."""
if rng is None:
rng = np.random.default_rng()
return (rng.random(size=n) < self._p).astype(int)
class DiscreteSampler(DiscreteDistribution):
"""
General discrete distribution from values and probabilities.
Allows specifying an arbitrary discrete distribution by providing
the possible values and their probabilities.
Supports two calling conventions (matching the MATLAB/JAR API):
- DiscreteSampler(p): values are implicitly 1, 2, ..., n (MATLAB-style)
- DiscreteSampler(p, x): explicit probabilities and values, in this
order (MATLAB DiscreteSampler(p,x) / JAR DiscreteSampler(p,x))
Args:
p: Array of probabilities (weights) of each item.
x: Array of possible values (only in the two-argument form).
"""
def __init__(self, p: Union[list, np.ndarray],
x: Optional[Union[list, np.ndarray]] = None):
super().__init__()
self._name = 'DiscreteSampler'
if x is None:
# MATLAB-style: DiscreteSampler(p) - values are 1, 2, ..., n
self._probs = np.array(p, dtype=float)
n = len(self._probs)
self._values = np.arange(1, n + 1, dtype=float) # 1-indexed like MATLAB
else:
# Full form: DiscreteSampler(p, x) - MATLAB/JAR argument order
self._probs = np.array(p, dtype=float)
self._values = np.array(x, dtype=float)
if len(self._values) != len(self._probs):
raise ValueError("values and probs must have the same length")
if len(self._values) == 0:
raise ValueError("Must have at least one value")
if np.any(self._probs < 0):
raise ValueError("Probabilities must be non-negative")
# Normalize probabilities
total = self._probs.sum()
if total <= 0:
raise ValueError("Probabilities must sum to a positive value")
self._probs = self._probs / total
# Sort by values for CDF computation
sort_idx = np.argsort(self._values)
self._values = self._values[sort_idx]
self._probs = self._probs[sort_idx]
self._cdf = np.cumsum(self._probs)
@property
def values(self) -> np.ndarray:
"""Get the possible values."""
return self._values.copy()
@property
def probs(self) -> np.ndarray:
"""Get the probabilities."""
return self._probs.copy()
def getMean(self) -> float:
"""Get the mean."""
return float(np.sum(self._values * self._probs))
def getVar(self) -> float:
"""Get the variance."""
mean = self.getMean()
return float(np.sum(self._probs * (self._values - mean) ** 2))
def getSupport(self) -> Tuple[float, float]:
"""Get the support [min, max]."""
return (float(self._values[0]), float(self._values[-1]))
def evalPMF(self, x: float) -> float:
"""Evaluate the PMF at point x."""
idx = np.where(np.isclose(self._values, x))[0]
if len(idx) == 0:
return 0.0
return float(self._probs[idx[0]])
def evalCDF(self, x: float) -> float:
"""Evaluate the CDF at point x."""
if x < self._values[0]:
return 0.0
idx = np.searchsorted(self._values, x, side='right') - 1
return float(self._cdf[min(idx, len(self._cdf) - 1)])
def sample(self, n: int = 1, rng: Optional[np.random.Generator] = None) -> np.ndarray:
"""Generate random samples."""
if rng is None:
rng = np.random.default_rng()
indices = rng.choice(len(self._values), size=n, p=self._probs)
return self._values[indices]
class EmpiricalCdf(DiscreteDistribution):
"""
Empirical distribution from CDF data.
The law is piecewise linear through the points (x_i, F_i), with an atom of
mass F_1 at x_1, matching MATLAB, the JAR and C++ (see
_kb/07-cross-language-parity.md, EmpiricalCDF semantics).
Args:
values: Nondecreasing support points x_i.
cdf: Cumulative probabilities F_i in [0, 1], nondecreasing, ending at 1.
"""
def __init__(self, values: Union[list, np.ndarray],
cdf: Union[list, np.ndarray]):
super().__init__()
# see _kb/07-cross-language-parity.md (EmpiricalCdf naming) for rationale
self._name = 'EmpiricalCdf'
x = np.asarray(values, dtype=float).ravel()
F = np.asarray(cdf, dtype=float).ravel()
if len(x) != len(F):
raise LineInputError("EmpiricalCdf: values and cdf must have the same length")
if len(x) == 0:
raise LineInputError("EmpiricalCdf: at least one point is required")
# Sort by (F, x) and drop duplicate rows, as MATLAB unique(...,'rows') does.
rows = np.unique(np.column_stack([F, x]), axis=0)
F, x = rows[:, 0], rows[:, 1]
tol = 1e-8
if np.any(F < 0) or np.any(F > 1 + tol):
raise LineInputError("EmpiricalCdf: cdf values must lie in [0, 1]; "
"the call is EmpiricalCdf(values, cdf)")
if np.any(np.diff(x) < 0):
raise LineInputError("EmpiricalCdf: values must be nondecreasing in the cdf")
if abs(F[-1] - 1.0) > tol:
raise LineInputError("EmpiricalCdf: the last cdf value must be 1")
self._values = x
self._cdf = F
@property
def values(self) -> np.ndarray:
"""Get the values."""
return self._values.copy()
@property
def cdf_values(self) -> np.ndarray:
"""Get the CDF values."""
return self._cdf.copy()
def _raw_moment(self, k: int) -> float:
"""E[X^k]: the atom F_1 x_1^k plus the midpoint rule over each linear segment."""
x, F = self._values, self._cdf
mid = (x[1:] + x[:-1]) / 2.0
return float(F[0] * x[0] ** k + np.sum(np.diff(F) * mid ** k))
def getMean(self) -> float:
"""Get the mean."""
return self._raw_moment(1)
def getVar(self) -> float:
"""Get the variance."""
m1 = self._raw_moment(1)
return self._raw_moment(2) - m1 ** 2
def getMoments(self):
"""Return (m1, m2, m3, scv, skew)."""
m1, m2, m3 = self._raw_moment(1), self._raw_moment(2), self._raw_moment(3)
var = m2 - m1 ** 2
scv = var / m1 ** 2 if m1 != 0 else float('inf')
skew = (m3 - 3 * m1 * var - m1 ** 3) / var ** 1.5 if var > 0 else 0.0
return m1, m2, m3, scv, skew
def getSupport(self) -> Tuple[float, float]:
"""Get the support [min, max]."""
return (float(self._values[0]), float(self._values[-1]))
def evalPMF(self, x: float) -> float:
"""Probability mass at x: only the atom F_1 at x_1 carries mass."""
return float(self._cdf[0]) if np.isclose(x, self._values[0]) else 0.0
def evalCDF(self, x: float) -> float:
"""Evaluate the CDF at point x by linear interpolation."""
xs, F = self._values, self._cdf
if x < xs[0]:
return 0.0
if x >= xs[-1]:
return 1.0
i = int(np.searchsorted(xs, x, side='right')) - 1
if xs[i + 1] == xs[i]:
return float(F[i + 1])
return float(F[i] + (F[i + 1] - F[i]) * (x - xs[i]) / (xs[i + 1] - xs[i]))
def evalLST(self, s: float) -> float:
"""Laplace-Stieltjes transform, consistent with the midpoint moments."""
x, F = self._values, self._cdf
mid = (x[1:] + x[:-1]) / 2.0
return float(F[0] * np.exp(-s * x[0]) + np.sum(np.diff(F) * np.exp(-s * mid)))
def sample(self, n: int = 1, rng: Optional[np.random.Generator] = None) -> np.ndarray:
"""Generate random samples by inverting the piecewise-linear CDF."""
if rng is None:
rng = np.random.default_rng()
u = rng.random(size=n)
xs, F = self._values, self._cdf
out = np.empty(n)
i = np.searchsorted(F, u, side='left')
atom = i == 0
out[atom] = xs[0]
j = np.minimum(i[~atom], len(F) - 1)
dF = F[j] - F[j - 1]
w = np.where(dF > 0, (u[~atom] - F[j - 1]) / np.where(dF > 0, dF, 1.0), 1.0)
out[~atom] = xs[j - 1] + w * (xs[j] - xs[j - 1])
return out
# see _kb/07-cross-language-parity.md (EmpiricalCdf naming) for rationale
EmpiricalCDF = EmpiricalCdf
[docs]
class Replayer(DiscreteDistribution):
"""
Trace-based distribution that replays recorded values.
Used for simulation where inter-arrival or service times are
read from a trace file or data array.
Args:
trace: Array of values to replay.
loop: Whether to loop when trace is exhausted (default: True).
"""
def __init__(self, trace: Union[str, list, np.ndarray], loop: bool = True):
super().__init__()
self._name = 'Replayer'
self._file_path = None # Store file path for JMT if provided
# Handle file path string
if isinstance(trace, str):
# Store file path for JMT
self._file_path = trace
# Read trace from file
self._trace = np.loadtxt(trace, dtype=float)
else:
self._trace = np.atleast_1d(np.array(trace, dtype=float))
if len(self._trace) == 0:
raise ValueError("Trace must have at least one value")
# see _kb/09-ldes-and-cache.md (Replayer/Trace, distributions save) for rationale
self._file_path_is_temp = False
if self._file_path is None:
# see _kb/09-ldes-and-cache.md (Replayer/Trace, distributions save) for rationale
import os
import tempfile
_fd, _tmp = tempfile.mkstemp(suffix='.trace', prefix='line_replayer_')
with os.fdopen(_fd, 'w') as _fh:
for _v in self._trace:
_fh.write('%.17g\n' % float(_v))
self._file_path = _tmp
self._file_path_is_temp = True
self._loop = loop
self._index = 0
@property
def trace(self) -> np.ndarray:
"""Get the trace data."""
return self._trace.copy()
@property
def loop(self) -> bool:
"""Check if looping is enabled."""
return self._loop
[docs]
def reset(self) -> None:
"""Reset the replay index to the beginning."""
self._index = 0
[docs]
def getMean(self) -> float:
"""Get the mean of the trace."""
return float(np.mean(self._trace))
[docs]
def getVar(self) -> float:
"""Get the variance of the trace."""
return float(np.var(self._trace))
[docs]
def getSkewness(self) -> float:
"""Get the skewness of the trace."""
from scipy.stats import skew
# Use bias=False for sample skewness (matches MATLAB's skewness(data,0))
return float(skew(self._trace, bias=False))
[docs]
def getSupport(self) -> Tuple[float, float]:
"""Get the support [min, max] of trace values."""
return (float(np.min(self._trace)), float(np.max(self._trace)))
[docs]
def evalPMF(self, x: float) -> float:
"""Evaluate PMF based on trace frequency."""
count = np.sum(np.isclose(self._trace, x))
return count / len(self._trace)
[docs]
def evalCDF(self, x: float) -> float:
"""Evaluate CDF based on trace."""
return np.sum(self._trace <= x) / len(self._trace)
[docs]
def evalLST(self, s: float) -> float:
"""Empirical Laplace-Stieltjes transform, mean(exp(-s*trace)).
Exact for the trace, and what the G/M/1 sigma-root needs: without it
sn.lst stayed None for a Replayer arrival and the MVA gm1 branch fell
back to a two-moment approximation. Matches MATLAB Replayer.evalLST and
the JAR Replayer.evalLST.
"""
return float(np.mean(np.exp(-s * self._trace)))
[docs]
def next_value(self) -> float:
"""Get the next value in the trace."""
if self._index >= len(self._trace):
if self._loop:
self._index = 0
else:
raise StopIteration("Trace exhausted and loop is disabled")
value = self._trace[self._index]
self._index += 1
return float(value)
[docs]
def sample(self, n: int = 1, rng: Optional[np.random.Generator] = None) -> np.ndarray:
"""Get the next n values from the trace."""
samples = np.zeros(n)
for i in range(n):
samples[i] = self.next_value()
return samples
[docs]
def isNHPP(self, **kwargs):
"""
Test whether the trace is a sample path of a NON-HOMOGENEOUS POISSON
process, by the conditional-uniform KS test with the Lewis refinement
(:func:`infer_nhpp_ks`).
WHY THE QUESTION IS WORTH ASKING. A Replayer is used wherever a measured
stream is fed to a solver, and every analytical method that consumes it
as an arrival process assumes SOMETHING about its dependence structure.
This test says whether the Poisson assumption -- independent increments,
whatever the rate does with time -- survives contact with the data, which
is the assumption a time-varying analysis (SolverFLD's ``mtginf``,
``mol``, ``tvms``) rests on. A small p-value says the stream is not
Poisson at any rate function, so those methods are answering a different
process.
The trace holds INTER-ARRIVAL times, so the arrival epochs are their
cumulative sum and the horizon is the last of them.
Returns:
The dict of :func:`infer_nhpp_ks`: ``statistic``, ``pvalue``, ``n``,
``uniforms`` and ``transformed``.
References:
S.-H. Kim, W. Whitt (2014). Are call center and hospital arrivals
well modeled by nonhomogeneous Poisson processes? Manufacturing &
Service Operations Management 16(3), 464-480.
"""
from ..api.infer import infer_nhpp_ks
ia = np.asarray(self.trace, dtype=float).ravel()
if ia.size < 2:
raise LineInputError('the trace needs at least two inter-arrival times to test', error_id='LINE:input:Argument')
epochs = np.cumsum(ia)
return infer_nhpp_ks(epochs, float(epochs[-1]), **kwargs)
is_nhpp = isNHPP
[docs]
def fit_exp(self):
"""
Fit an exponential to the trace mean (MATLAB Replayer.fitExp).
Returns:
Exp distribution with the mean of the trace.
"""
from .continuous import Exp
return Exp.fit_mean(self.getMean())
[docs]
def fit_coxian(self):
"""
Fit a two-phase Coxian to the first three moments of the trace
(MATLAB Replayer.fitCoxian, which calls Cox2.fit).
Returns:
Cox2 distribution fitted to the trace moments.
"""
from .markovian import Cox2
mean = self.getMean()
var = self.getVar()
return Cox2.fit_central(mean, var, self.getSkewness())
[docs]
def fitExp(self):
"""camelCase alias of fit_exp (MATLAB/JAR spelling)."""
return self.fit_exp()
[docs]
def fitCoxian(self):
"""camelCase alias of fit_coxian (MATLAB/JAR spelling)."""
return self.fit_coxian()
[docs]
def fitAPH(self):
"""camelCase alias of fit_aph (MATLAB/JAR spelling)."""
return self.fit_aph()
[docs]
def fit_aph(self):
"""
Fit an acyclic phase-type (APH) distribution to the trace data.
Uses 3-moment matching (mean, SCV, skewness) to determine the
optimal APH representation, matching MATLAB's behavior.
Returns:
APH distribution fitted to the trace data.
"""
from .markovian import APH
mean = self.getMean()
var = self.getVar()
scv = var / (mean ** 2) if mean > 0 else 1.0
skewness = self.getSkewness()
# Compute raw moments from central moments
e1 = mean
e2 = (1 + scv) * e1 ** 2
e3 = -(2 * e1 ** 3 - 3 * e1 * e2 - skewness * (e2 - e1 ** 2) ** (3 / 2))
# Use APHFrom3Moments for 3-moment matching (matches MATLAB)
try:
alpha, T = self._aph_from_3_moments([e1, e2, e3])
# Convert from matrix to array format
alpha = np.asarray(alpha).flatten()
T = np.asarray(T)
return APH(alpha, T)
except (ValueError, ArithmeticError):
# Fall back to 2-moment matching if 3-moment fails
return APH(mean=mean, scv=scv)
@staticmethod
def _aph_from_3_moments(moms, maxSize=100):
"""
Returns an acyclic PH which has the same 3 moments as given.
Determines the order and structure automatically to match the given moments.
Based on BUTools APHFrom3Moments.
"""
import numpy.matlib as ml
import math
import cmath
def _aph_2nd_moment_lower_bound(m1, n):
return float(m1) * m1 * (n + 1) / n
def _aph_3rd_moment_lower_bound(m1, m2, n):
n2 = m2 / m1 / m1
if n2 < (n + 1.0) / n:
return np.inf
elif n2 < (n + 4.0) / (n + 1.0):
p = ((n + 1.0) * (n2 - 2.0)) / (3.0 * n2 * (n - 1.0)) * (
(-2.0 * math.sqrt(n + 1.0)) / cmath.sqrt(-3.0 * n * n2 + 4.0 * n + 4.0) - 1.0)
a = (n2 - 2.0) / (p * (1.0 - n2) + cmath.sqrt(p * p + p * n * (n2 - 2.0) / (n - 1.0)))
l = ((3.0 + a) * (n - 1.0) + 2.0 * a) / ((n - 1.0) * (1.0 + a * p)) - (2.0 * a * (n + 1.0)) / (
2.0 * (n - 1.0) + a * p * (n * a + 2.0 * n - 2.0))
return l.real * m1 * m2
else:
return (n + 1.0) / n * n2 * m1 * m2
def _aph_3rd_moment_upper_bound(m1, m2, n):
n2 = m2 / m1 / m1
if n2 < (n + 1.0) / n:
return -np.inf
elif n2 <= n / (n - 1.0):
return m1 * m2 * (2.0 * (n - 2.0) * (n * n2 - n - 1.0) * math.sqrt(1.0 + (n * (n2 - 2.0)) / (n - 1.0)) + (
n + 2.0) * (3.0 * n * n2 - 2.0 * n - 2.0)) / (n * n * n2)
else:
return np.inf
def _norm_moms_from_moms(m):
return [float(m[i]) / m[i - 1] / m[0] if i > 0 else m[0] for i in range(len(m))]
m1, m2, m3 = moms
# detect number of phases needed
n = 2
while n < maxSize and (_aph_2nd_moment_lower_bound(m1, n) > m2 or _aph_3rd_moment_lower_bound(m1, m2, n) >= m3 or _aph_3rd_moment_upper_bound(m1, m2, n) <= m3):
n = n + 1
# if PH is too large, adjust moment to bounds
if _aph_2nd_moment_lower_bound(m1, n) > m2:
m2 = _aph_2nd_moment_lower_bound(m1, n)
if _aph_3rd_moment_lower_bound(m1, m2, n) > m3:
m3 = _aph_3rd_moment_lower_bound(m1, m2, n)
if _aph_3rd_moment_upper_bound(m1, m2, n) < m3:
m3 = _aph_3rd_moment_upper_bound(m1, m2, n)
# compute normalized moments
n1, n2, n3 = _norm_moms_from_moms([m1, m2, m3])
if n2 > 2.0 or n3 < 2.0 * n2 - 1.0:
b = (2.0 * (4.0 - n * (3.0 * n2 - 4.0)) / (n2 * (4.0 + n - n * n3) + math.sqrt(n * n2) * math.sqrt(
12.0 * n2 * n2 * (n + 1.0) + 16.0 * n3 * (n + 1.0) + n2 * (n * (n3 - 15.0) * (n3 + 1.0) - 8.0 * (n3 + 3.0))))).real
a = (b * n2 - 2.0) * (n - 1.0) * b / (b - 1.0) / n
p = (b - 1.0) / a
lamb = (p * a + 1.0) / n1
mu = (n - 1.0) * lamb / a
# construct representation
alpha = ml.zeros((1, n))
alpha[0, 0] = p
alpha[0, n - 1] = 1.0 - p
A = ml.zeros((n, n))
A[n - 1, n - 1] = -lamb
for i in range(n - 1):
A[i, i] = -mu
A[i, i + 1] = mu
return (alpha, A)
else:
c4 = n2 * (3.0 * n2 - 2.0 * n3) * (n - 1.0) * (n - 1.0)
c3 = 2.0 * n2 * (n3 - 3.0) * (n - 1.0) * (n - 1.0)
c2 = 6.0 * (n - 1.0) * (n - n2)
c1 = 4.0 * n * (2.0 - n)
c0 = n * (n - 2.0)
fs = np.roots([c4, c3, c2, c1, c0])
for f in fs:
if abs((n - 1) * (n2 * f * f * 2 - 2 * f + 2) - n) < 1e-14:
continue
a = 2.0 * (f - 1.0) * (n - 1.0) / ((n - 1.0) * (n2 * f * f - 2.0 * f + 2.0) - n)
p = (f - 1.0) * a
lamb = (a + p) / n1
mu = (n - 1.0) / (n1 - p / lamb)
if np.isreal(p) and np.isreal(lamb) and np.isreal(mu) and p >= 0 and p <= 1 and lamb > 0 and mu > 0:
alpha = ml.zeros((1, n))
alpha[0, 0] = p.real
alpha[0, 1] = 1.0 - p.real
A = ml.zeros((n, n))
A[0, 0] = -lamb.real
A[0, 1] = lamb.real
for i in range(1, n):
A[i, i] = -mu.real
if i < n - 1:
A[i, i + 1] = mu.real
return (alpha, A)
raise ValueError("No APH found for the given 3 moments!")
class Trace(Replayer):
"""
Empirical time series from a trace file.
Alias for Replayer with additional moment computation for
histogram-style trace data (value, count pairs).
Args:
data: Array of values to replay, or 2-column array of (x, cdf) pairs.
loop: Whether to loop when trace is exhausted (default: True).
"""
def __init__(self, data: Union[list, np.ndarray], loop: bool = True):
data = np.atleast_2d(np.array(data, dtype=float))
# If 2-column data (x, cdf format), extract values
if data.shape[1] == 2:
self._histogram_data = data.copy()
# Convert CDF to samples (approximate)
values = data[:, 0]
super().__init__(values, loop)
else:
self._histogram_data = None
super().__init__(data.flatten(), loop)
self._name = 'Trace'
def getMoments(self) -> Tuple[float, float, float, float, float]:
"""
Compute moments from histogram-style trace data.
Returns:
Tuple of (m1, m2, m3, scv, skew) - first three moments, SCV, and skewness.
"""
if self._histogram_data is not None:
data = self._histogram_data
m1 = 0.0 # First moment
m2 = 0.0 # Second moment
m3 = 0.0 # Third moment
for i in range(len(data) - 1):
# Trapezoidal integration
mid_val = (data[i + 1, 1] - data[i, 1]) / 2 + data[i, 1]
width = data[i + 1, 0] - data[i, 0]
m1 += mid_val * width
m2 += mid_val ** 2 * width
m3 += mid_val ** 3 * width
if m1 > 0:
scv = (m2 / m1 ** 2) - 1
var = m2 - m1 ** 2
if var > 0:
skew = (m3 - 3 * m1 * var - m1 ** 3) / (var ** 1.5)
else:
skew = 0.0
else:
scv = 0.0
skew = 0.0
return (m1, m2, m3, scv, skew)
else:
# Compute from raw trace data
m1 = float(np.mean(self._trace))
m2 = float(np.mean(self._trace ** 2))
m3 = float(np.mean(self._trace ** 3))
var = m2 - m1 ** 2
scv = var / m1 ** 2 if m1 > 0 else 0.0
skew = (m3 - 3 * m1 * var - m1 ** 3) / (var ** 1.5) if var > 0 else 0.0
return (m1, m2, m3, scv, skew)