me_sample.h File Reference

Sample a matrix exponential by numerical inversion of its exact CDF. More...

#include <cmath>
#include <cstddef>
#include <vector>
#include "line/api/mam/map_moment.h"
#include "line/api/pfqn/pfqn_mc_common.h"
#include "line/num/number.h"
#include "line/util/eig.h"
#include "line/util/error.h"
#include "line/util/expm.h"
#include "line/util/matrix.h"

Include dependency graph for me_sample.h:

Go to the source code of this file.

Classes
class	line::mam::MeSampler< T >
	Stateful ME sampler holding the inversion table. More...

Namespaces
namespace	line
namespace	line::mam

Functions
template<class T>
std::vector< T >	line::mam::me_sample (const Map< T > &m, std::size_t n, pfqn::McRng &rng, const std::vector< T > &a0=std::vector< T >())
	me_sample: n INDEPENDENT variates of the matrix exponential (D0, D1).

Detailed Description

Sample a matrix exponential by numerical inversion of its exact CDF.

Templated port of matlab/lib/kpctoolbox/map/me_sample.m, cross-checked against jar/src/main/java/jline/api/mam/Me_sample.java.

AN ME IS A RENEWAL PROCESS, which is what separates this from rap_sample. The inter-arrival times are i.i.d. draws from F(t) = 1 - pie exp(D0 t) e, so the entry law is the SAME at every sample and the inversion table built from it can be reused. rap_sample instead advances the entry law after each arrival, which is correct for a RAP and wrong for an ME: it delivers a correlated stream where the reference delivers an independent one.

WHY A TABLE AND NOT A ROOT SOLVE. Routing ME through rap_sample costs a full matrix exponential per bisection step, some fifty of them per variate; on a 15-phase CME that is 140x the Java engine on the same model, enough for the LDES wrapper to hit its wall-clock budget and return an empty result. The reference tabulates the CDF once on a grid whose horizon is doubled until the survival is below TAILTOL, locates the variate by binary search and polishes it with Newton steps on the exact CDF and density, so the answer is not limited by the linear interpolation of the table.

ONE UNIFORM PER VARIATE, as in rap_sample, so a caller switching between the two consumes the generator at the same rate.

The grid itself is generated by a single expm(D0 h) and repeated vector-matrix products, and the Newton step propagates the row vector by a scaled Taylor series rather than re-exponentiating: both are O(K^2) where a matrix exponential is O(K^3).

ARITHMETIC: transcendental.

Definition in file me_sample.h.