ctmc_fau.h File Reference

Transient distribution of a CTMC by fast adaptive uniformization. More...

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <iterator>
#include <vector>
#include "line/api/mc/ctmc_foxglynn.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for ctmc_fau.h:

Go to the source code of this file.

Classes
struct	line::mc::FauResult< T >

Namespaces
namespace	line
namespace	line::mc

Functions
template<class T>
FauResult< T >	line::mc::ctmc_fau (const std::vector< T > &pi0, const Matrix< T > &Q, const T &t, double epsilon=1e-6, double delta=1e-12, long maxsteps=-1)
	Transient distribution of a CTMC by fast adaptive uniformization.

Variables
constexpr long	line::mc::FAU_MAX_STEPS = 1000000
	Default cap on birth steps, so a pathological horizon reports truncation.

Detailed Description

Transient distribution of a CTMC by fast adaptive uniformization.

Templated port of matlab/src/api/mc/ctmc_fau.m and jar/src/main/java/jline/api/mc/Ctmc_fau.java.

Ordinary uniformization fixes one rate q >= max_i |q_ii| over the WHOLE state space and mixes the powers of P = I + Q/q against a Poisson(q t) law, so its cost is set by the fastest state anywhere, including states carrying no probability at time t. Adaptive uniformization (van Moorsel and Sanders, 1994) picks a rate per step from the states the iterate occupies,

Lambda_n >= max{|q_ii| : i in supp(u^(n))}, u^(n+1) = u^(n)(I + Q/Lambda_n),

which keeps every entry of u^(n+1) nonnegative. The subordinating process is then the pure birth process N(t) with rates Lambda_0, Lambda_1, ... and pi(t) = sum_n P{N(t) = n} u^(n). The fast variant (Mateescu, Wolf, Didier and Henzinger, 2010) drops an entry below delta rather than propagating it, so the support tracks the states of non-negligible occupancy instead of the reachable set.

NOTHING IS RENORMALIZED, so the error is measured rather than estimated: the birth index truncated at K, the Poisson window of the weight computation and the delta threshold each remove mass and none puts any back, whence 0 <= pi(t) - pit componentwise and |pi(t) - pit|_1 = sum(pi0) - sum(pit), which is what errorBound reports.

The birth weights are exact rather than quadratured: the rates generate a bidiagonal generator on the birth index plus one absorbing overflow index, and its transient distribution is obtained by uniformizing that scalar chain at Lstar = max_n Lambda_n and mixing the shipped Fox-Glynn weights. The sweep runs twice because b_n(t) needs the rates up to n, which are not known before the sweep ends, while u^(n) is needed after them, and storing every iterate would cost K times the support.

GATED ON TRANSCENDENTAL ARITHMETIC, exactly as ctmc_foxglynn is: the stopping rule is a logarithmic tail bound and the mixture is an approximation of exp(Qt) controlled by epsilon, so the rate bookkeeping and the truncation decisions are taken in double while the iterate and the weighted sum run in T, which is what a high-precision instantiation buys.

This is a transient method: it produces no stationary distribution.

Definition in file ctmc_fau.h.