ctmc_saddlepoint.h File Reference

Saddlepoint approximation of Pr{N(t)=k} for the counting process of a MAP. More...

#include <algorithm>
#include <cmath>
#include <complex>
#include <limits>
#include <cstddef>
#include <numeric>
#include <string>
#include <vector>
#include "line/num/number.h"
#include "line/util/eig.h"
#include "line/util/error.h"
#include "line/util/expm.h"
#include "line/api/mc/ctmc_solve.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"

Include dependency graph for ctmc_saddlepoint.h:

Go to the source code of this file.

Classes
struct	line::mc::PerronState< T >
	Perron root of A(theta) with its first two derivatives and the amplitude. More...
struct	line::mc::SaddlepointResult< T >
	One entry per (t,k) pair. More...

Namespaces
namespace	line
namespace	line::mc

Enumerations
enum	line::mc::SaddlepointMethod { line::mc::SADDLEPOINT_DANIELS2 = 0 , line::mc::SADDLEPOINT_DANIELS = 1 , line::mc::SADDLEPOINT_PLAIN = 2 }
	Which term of the steepest-descent expansion to stop at. More...

Functions
template<class T>
PerronState< T >	line::mc::ctmc_saddlepoint_perron (const Matrix< T > &D0, const Matrix< T > &D1, const std::vector< T > &pi0, const T &th)
	Perron root of A(th) = D0 + exp(th)*D1 with deta, d2eta and the amplitude.
template<class T>
SaddlepointResult< T >	line::mc::ctmc_saddlepoint (const Matrix< T > &D0, const Matrix< T > &D1, const std::vector< T > &t, const std::vector< long > &k, SaddlepointMethod method=SADDLEPOINT_DANIELS2, const std::vector< T > &pi0=std::vector< T >())
	Pr{N(t)=k} over arrays of horizons and counts.
template<class T>
T	line::mc::ctmc_saddlepoint (const Matrix< T > &D0, const Matrix< T > &D1, const T &t, long k, SaddlepointMethod method=SADDLEPOINT_DANIELS2)
	Pr{N(t)=k} at a single (t,k), with the default method.
template<class T>
T	line::mc::tx_exp (const T &x)
	exp/log/isfinite through ADL, the idiom the rest of api/mc uses: using std::exp then an unqualified call, so a Real<D> or Rational picks up its own overload.
template<class T>
T	line::mc::tx_log (const T &x)
template<class T>
bool	line::mc::tx_finite (const T &x)

Variables
static const double	line::mc::CTMC_SADDLEPOINT_K2_MIN = 5.0
	Below this value of K2 = t*eta''(theta*) the expansion is out of its regime.

Detailed Description

Saddlepoint approximation of Pr{N(t)=k} for the counting process of a MAP.

Templated port of matlab/src/api/mc/ctmc_saddlepoint.m. 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 result carries a flag saying so.

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.

GATED ON TRANSCENDENTAL ARITHMETIC (exp, log, sqrt) and, for the Perron root, on LAPACK through eig_values, as several other api/mam headers already are.

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.

Definition in file ctmc_saddlepoint.h.