Package jline.api.mc

Class Ctmc_saddlepoint

java.lang.Object
jline.api.mc.Ctmc_saddlepoint
public final class Ctmc_saddlepoint extends Object
Saddlepoint approximation of Pr{N(t)=k} for the counting process of a MAP.

The probability that the Markovian arrival process (D0,D1) records exactly k events in (0,t], obtained by steepest-descent inversion of the counting generating function instead of by forming the k-th superdiagonal block of expm(t*X).

The counting generating function is the matrix exponential 

     sum_k P(k,t) z^k = expm(t*(D0 + z*D1)),
so the cumulant generating function of N(t) is eta(theta) = spectral abscissa of A(theta) = D0 + exp(theta)*D1, the Perron root of an irreducible Metzler matrix: real, simple, strictly convex in theta, with eta(0)=0 and eta'(0)=lambda. Inverting by steepest descent gives Daniels (1954), 
     Pr{N(t)=k} ~ g(theta*) * exp(t*eta(theta*) - k*theta*)
                            / sqrt(2*pi*t*eta''(theta*)),
with the saddle theta* solving eta'(theta*) = k/t and g the amplitude of the Perron projection, g(theta) = (pi0*v)*(u*1), u and v the left and right Perron vectors normalised by u*v = 1.

THE EXPANSION PARAMETER IS K2 = t*eta''(theta*), THE VARIANCE OF THE COUNT, not its mean and not t. Measured error laws, with the constants flat to two digits over Erlang orders 1..8 and horizons 10..160: 

     err(DANIELS) = 0.083 / K2        err(DANIELS2) = 0.017 / K2^2
For a renewal Erlang(r) the count variance rate is lambda/r, so K2 = lambda*t/r and an Erlang-4 at t=50 is as accurate as a Poisson at t=12.5: low variability shrinks the parameter, it does not break the method. Below K2 = 5 the expansion is out of its regime and the call warns.

This is an asymptotic method, not a quadrature: use it for rare-event and large-deviation coefficients, where k/t is away from lambda or where the probability underflows. For the bulk of the transient distribution, i.e. every block k=0..N-1 at once at moderate t, uniformization (Ctmc_uniformization, Ctmc_foxglynn) is both exact and faster.

ATTRIBUTION. The first-order form is Daniels (1954). The amplitude g and the whole 'daniels2' bracket are NOT a rederivation: they are Jensen, "Saddlepoint Expansions for Sums of Markov Dependent Variables on a Continuous State Space", Probab. Th. Rel. Fields 89, 1991, Eq. (4.4) with the coefficients on p.191. His gamma_0(s) = (sum_i c_i)(sum_i r_i P(Y_0=i)) is exactly g under his own normalisation sum_i r_i c_i = 1, and expanding his alpha_0 + (1/n){-alpha_3/2 + alpha_4/8 - 5*alpha_5/24} reproduces g*(1 + lam4/8 - 5*lam3^2/24) - g''/(2*K2) + g'*K3/(2*K2^2) term for term; his Theorem 4.1 gives the O(n^-2) error measured here as 0.017/K2^2. Jensen works with discrete-n sums over a Markov chain, so the continuous-time MAP counting process is that result transcribed, n -> t and the kernel eigenvalue -> the Perron root of D0+exp(theta)*D1.

Since:
LINE 3.0

  • Field Details

    • K2_MIN

      public static final double K2_MIN
      Below this value of K2 = t*eta''(theta*) the expansion is out of its regime. Do NOT threshold on lambda*t: for Erlang(r) the count variance rate is lambda/r, so K2 = lambda*t/r, and lambda*t over-warns on Poisson-like processes while under-warning on low-variability ones.
      See Also:

    • DANIELS2

      public static final String DANIELS2
      Second-order saddlepoint, error O(1/K2^2). The default.
      See Also:

    • DANIELS

      public static final String DANIELS
      First-order saddlepoint with the Perron amplitude, error O(1/K2).
      See Also:

    • PLAIN

      public static final String PLAIN
      Bare first-order form with the amplitude set to 1.
      See Also:

  • Method Details

    • ctmc_saddlepoint

      public static double ctmc_saddlepoint(Matrix D0, Matrix D1, double t, int k)
      Pr{N(t)=k} at a single (t,k), with the default method.
      Parameters:
      D0 - generator of the phase process with the counted transitions removed
      D1 - rates of the counted transitions; D0+D1 must be an irreducible generator
      t - time horizon
      k - event count, a nonnegative integer
      Returns:
      the approximation

    • ctmc_saddlepoint

      public static double ctmc_saddlepoint(Matrix D0, Matrix D1, double t, int k, String method)
      Pr{N(t)=k} at a single (t,k).
      Parameters:
      D0 - generator of the phase process with the counted transitions removed
      D1 - rates of the counted transitions
      t - time horizon
      k - event count
      method - DANIELS2 (default), DANIELS or PLAIN
      Returns:
      the approximation

    • ctmc_saddlepoint

      public static Ctmc_saddlepoint.SaddlepointResult ctmc_saddlepoint(Matrix D0, Matrix D1, double[] t, int[] k, String method, double[] pi0)
      Pr{N(t)=k} over arrays of horizons and counts.
      Parameters:
      D0 - generator of the phase process with the counted transitions removed
      D1 - rates of the counted transitions
      t - time horizons; length 1 broadcasts against k
      k - event counts; length 1 broadcasts against t
      method - DANIELS2 (default), DANIELS or PLAIN; null selects the default
      pi0 - initial phase distribution; null selects the stationary distribution of D0+D1
      Returns:
      the per-point approximations and diagnostics