Sample path of a continuous-time Markov chain given its generator. More...

#include <cmath>
#include <cstddef>
#include <vector>
#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 ctmc_simulate.h:

Classes
struct	line::mc::CtmcPath< T >
	One simulated sample path: the state visited at each step and its holding time. More...

Namespaces
namespace	line
namespace	line::mc

Functions
template<class T>
std::vector< T >	line::mc::ctmc_simulate_reference_initial_law (const std::vector< T > &pi0)
	The initial-state law that ctmc_simulate.m actually realizes for a given pi0, as opposed to pi0 itself.
template<class T>
CtmcPath< T >	line::mc::ctmc_simulate (const Matrix< T > &Q, const std::vector< T > &pi0, std::size_t n, pfqn::McRng &rng)
	Simulate n steps of the CTMC with generator Q.

Detailed Description

Sample path of a continuous-time Markov chain given its generator.

Templated port of matlab/src/api/mc/ctmc_simulate.m. The chain is simulated by the standard jump-chain construction: from state i the holding time is exponential with mean -1/Q(i,i), and the next state is drawn from the embedded jump chain P(i,j) = Q(i,j) / sum_{k != i} Q(i,k).

This is the one entry point in this batch that takes MATRICES and not an sn: its arguments are the generator, an initial distribution and a step count, so it is portable independently of the NetworkStruct layer.

ARITHMETIC. Gated on num_traits<T>::has_transcendental. The holding times are exponential deviates, drawn by inverse transform as -mean * log(u), so there is no exact instantiation: the sample path is a realization, not a number that a rational field could represent.

RANDOMNESS. The generator is line::pfqn::McRng (a std::mt19937_64) passed by reference and advanced by the call, the same convention every Monte Carlo entry point in this tree uses, so a caller controls reproducibility by controlling the generator state. The stream is NOT comparable with MATLAB's – different generator, different mapping from bits to deviates – so the oracle for this function is distributional, never path-for-path. Two generator steps are consumed per simulated step: one for the holding time and one for the jump.

REFERENCE DEFECT: the initial state is not drawn from pi0.

ctmc_simulate.m selects the starting state with

[~, st] = min(abs(rand - cumsum(pi0)));

which returns the state whose CUMULATIVE probability is nearest to the uniform deviate, ties going to the lowest index because min returns the first minimizer. That partitions [0,1] at the MIDPOINTS between consecutive DISTINCT cumulative values c_k = sum_{j<=k} p_j, instead of at the cumulative values themselves. Where the c_k are distinct this reduces to

P(1) = p_1 + p_2/2,   P(k) = (p_k + p_{k+1})/2,   P(n) = p_n/2,

i.e. the LAST state always receives about half its intended mass and the first receives an excess. Where two consecutive c_k coincide – which is exactly what a zero-probability entry produces – the first of them takes the whole window and the rest get nothing, so the naive midpoint reading above does NOT apply there. Measured in MATLAB over 400000 draws:

pi0 = [0.5 0.5] -> [0.7498 0.2502] (want [0.5 0.5]) pi0 = [1/3 1/3 1/3] -> [0.4997 0.3335 0.1668] (want [1/3 1/3 1/3]) pi0 = [0.9 0 0.1] -> [0.9497 0.0000 0.0503] (want [0.9 0 0.1])

the last of which is also reproducible through the entry point itself, ctmc_simulate(Q, [0.9;0;0.1], 1) over 20000 calls giving [0.9488 0.0000 0.0512]. The error washes out of a long ergodic run – the time-average occupancy of a two-state chain still converges to the exact stationary law – so it is invisible in steady-state use and corrupts exactly the transient and short-run use that passing pi0 is for.

The port draws the initial state by correct inverse transform. Reproducing the defect was rejected: a sampler that does not sample from the distribution it is handed has no contract left to preserve, and unlike a closed-form value there is nothing downstream that could be calibrated against the wrong answer. ctmc_simulate_reference_initial_law below returns the law the reference actually realizes, so the discrepancy is available to a caller as a value rather than only as prose.

A second, smaller divergence: an ABSORBING state (a row whose off-diagonal entries are all zero) gives MATLAB F = 0/0 = NaN on that row, after which find(rand - NaN > 0) is empty and the chain silently jumps to state 1, while the holding time is exprnd(Inf) = Inf. The port raises NumericError naming the state instead: there is no Inf in an exact field, and a silent jump to state 1 is not a property of the chain.

Definition in file ctmc_simulate.h.

Classes

Namespaces

Functions

Detailed Description