map_sample.h File Reference

Sample the inter-arrival times of a MAP, a RAP or a matrix exponential. 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/error.h"
#include "line/util/matrix.h"

Include dependency graph for map_sample.h:

Go to the source code of this file.

Classes
struct	line::mam::SampleTrace
	The state a sampled inter-arrival began and ended in. More...

Namespaces
namespace	line
namespace	line::mam

Functions
template<class T>
std::size_t	line::mam::randp (const std::vector< T > &P, pfqn::McRng &rng)
	Draw an index from a discrete law, m3a's randp.
template<class T>
std::vector< T >	line::mam::map_sample (const Map< T > &m, std::size_t n, pfqn::McRng &rng, const std::vector< T > &pie0=std::vector< T >(), SampleTrace *trace=0)
	Sample the inter-arrival times of a MAP, a RAP or a matrix exponential.
template<class T>
std::vector< T >	line::mam::rap_sample (const Map< T > &m, std::size_t n, pfqn::McRng &rng, const std::vector< T > &a0=std::vector< T >(), std::vector< T > *a_out=0)
	Sample a RAP or a matrix exponential by inverse transform.

Detailed Description

Sample the inter-arrival times of a MAP, a RAP or a matrix exponential.

Templated port of matlab/lib/kpctoolbox/map/map_sample.m, rap_sample.m and me_sample.m, together with m3a's randp.m.

A MAP is simulated on its own state space: from the current phase the process takes hidden transitions (D0's off-diagonal) until one of the arrival transitions (D1) fires, and the inter-arrival time is the sum of the exponential holding times spent on the way. The reference draws the whole path and then adds the holding times; this port adds them as it goes, which is the same sum with no path buffer and no growing visits matrix.

A RAP IS NOT SIMULATED THAT WAY. Its D0 has negative off-diagonals, so there is no embedded jump chain to walk: the "phase" is a signed vector, not a state. The reference samples it by INVERSE TRANSFORM on the conditional distribution, advancing the row vector a <- a exp(D0 x) D1 / (a exp(D0 x) D1 e) after each arrival at x. That is what rap_sample does here, and it is why it costs a matrix exponential per sample where map_sample costs a handful of deviates. An ME is a RENEWAL process, so it does not pay that price: it lives in me_sample.h, which tabulates its CDF once.

RANDOMNESS. The generator is line::pfqn::McRng passed by reference and advanced by the call, the convention every Monte Carlo entry point in this tree uses. The stream is NOT comparable with MATLAB's – different generator, different mapping from bits to deviates – so the oracle for these functions is distributional, never path-for-path.

ARITHMETIC: transcendental. Exponential deviates and, for the RAP and ME paths, a matrix exponential.

Definition in file map_sample.h.