Class ME

All Implemented Interfaces:
Serializable, Copyable
Direct Known Subclasses:
CME

public class ME extends Markovian
A 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). Representation: - alpha: initial vector (may have negative entries or sum != 1) - A: matrix parameter (must have all eigenvalues with negative real parts) - Dominant eigenvalue of A must be negative and real The distribution function is: F(t) = 1 - alpha * exp(A*t) * e (under certain conditions) Moments: m_k = k! * alpha * (-A)^(-k) * e
See Also:
  • Constructor Details

    • ME

      public ME(Matrix alpha, Matrix A)
      Creates a Matrix Exponential distribution with specified initial vector and matrix parameter.
      Parameters:
      alpha - the initial vector (row vector as Matrix)
      A - the matrix parameter (must be square with negative real eigenvalues)
      Throws:
      IllegalArgumentException - if the representation is invalid
    • ME

      protected ME(Matrix alpha, Matrix A, boolean checkDensity)
      Creates a Matrix Exponential distribution, optionally skipping the density scan.
      Parameters:
      alpha - the initial vector (row vector as Matrix)
      A - the matrix parameter (must be square with negative real eigenvalues)
      checkDensity - scan the density for a negative value. 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.
      Throws:
      IllegalArgumentException - if the representation is invalid
  • Method Details

    • scanNegativeDensity

      public static ME.NegativeDensityScan scanNegativeDensity(Matrix alphaRow, Matrix A)
      Searches 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 is used in place of 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 ME_SCAN_TAIL 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.

      Parameters:
      alphaRow - the initial vector, as a 1 x n row
      A - the matrix parameter
      Returns:
      the scan outcome
    • getAlpha

      public Matrix getAlpha()
      Gets the initial vector alpha.
      Returns:
      the initial vector as a Matrix
    • getA

      public Matrix getA()
      Gets the matrix parameter A.
      Returns:
      the matrix parameter A
    • getNumberOfPhases

      public long getNumberOfPhases()
      Description copied from class: Markovian
      Gets the number of phases in this Markovian distribution.
      Overrides:
      getNumberOfPhases in class Markovian
      Returns:
      the number of phases
    • getProcess

      public MatrixCell getProcess()
      Description copied from class: Markovian
      Gets the matrix representation of this Markovian process.
      Overrides:
      getProcess in class Markovian
      Returns:
      MatrixCell containing D0, D1, ... matrices
    • sample

      public double[] sample(int n, Random random)
      Generates samples from this ME distribution.

      Sampling inverts the exact CDF F(t) = 1 - alpha*exp(A*t)*e. The CTMC walk inherited from Markovian (map_sample) is not used because it presumes a phase-type interpretation of (alpha, A), which fails whenever alpha has negative entries or A has negative off-diagonal entries.

      Overrides:
      sample in class Markovian
      Parameters:
      n - the number of samples to generate
      random - the random number generator to use
      Returns:
      array of n i.i.d. samples
    • getMean

      public double getMean()
      Gets the mean, m1 = -alpha*A^(-1)*e.

      Overrides the Markovian formula 1/map_lambda, which obtains the rate from the stationary vector of D0+D1. That vector is a probabilistic object of a Markovian process, and computing it for an ME means solving a linear system whose conditioning degrades with the oscillation of A: a CME of order 101 came out with a relative error of 4e-8, where the definition below is exact to 1e-13. Native Python ME.getMean already uses the definition.

      Overrides:
      getMean in class Markovian
      Returns:
      the mean
    • getVar

      public double getVar()
      Gets the variance, m2 - m1^2 with m2 = 2*alpha*A^(-2)*e.
      Overrides:
      getVar in class Markovian
      Returns:
      the variance
    • getSCV

      public double getSCV()
      Gets the squared coefficient of variation, var/mean^2.
      Overrides:
      getSCV in class Markovian
      Returns:
      the squared coefficient of variation
    • fitMoments

      public static ME fitMoments(double[] moments)
      Creates an ME distribution by fitting the given moments. Uses BuTools MEFromMoments algorithm.
      Parameters:
      moments - array of moments (requires 2*M-1 moments for order M ME distribution)
      Returns:
      an ME distribution matching the given moments
      Throws:
      IllegalArgumentException - if moments are invalid or fitting fails
    • fromExp

      public static ME fromExp(double rate)
      Creates an ME distribution from an exponential distribution. This is a convenience method showing that Exp is a special case of ME.
      Parameters:
      rate - the rate parameter (lambda)
      Returns:
      an ME distribution equivalent to Exp(rate)
    • fromErlang

      public static ME fromErlang(int k, double rate)
      Creates an ME distribution from an Erlang distribution. This is a convenience method showing that Erlang is a special case of ME.
      Parameters:
      k - number of phases
      rate - rate parameter for each phase
      Returns:
      an ME distribution equivalent to Erlang(k, rate)
    • fromHyperExp

      public static ME fromHyperExp(double[] p, double[] rates)
      Creates an ME distribution from a HyperExponential distribution. This is a convenience method showing that HyperExp is a special case of ME.
      Parameters:
      p - array of probabilities for each branch
      rates - array of rates for each branch
      Returns:
      an ME distribution equivalent to HyperExp(p, rates)