pfqn_mci.h File Reference

Monte Carlo Integration estimate of the normalizing constant of a closed product-form network (Ross, Wang and Yao; MonteQueue 2.0). More...

#include <cmath>
#include <cstddef>
#include <vector>
#include "line/api/pfqn/pfqn_bs.h"
#include "line/api/pfqn/pfqn_ca.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 pfqn_mci.h:

Go to the source code of this file.

Namespaces
namespace	line
namespace	line::pfqn

Enumerations
enum class	line::pfqn::MciVariant { line::pfqn::Imci , line::pfqn::Mci , line::pfqn::Rm }
	The proposal-rate rules the reference selects between. More...

Functions
template<class T>
NcResult< T >	line::pfqn::pfqn_mci (const Matrix< T > &D, const std::vector< int > &N, const std::vector< T > &Z, std::size_t samples, MciVariant variant, McRng &rng)
	Monte Carlo Integration estimate of the normalizing constant of a closed product-form network (Ross, Wang and Yao; MonteQueue 2.0).
template<class T>
NcResult< T >	line::pfqn::pfqn_mci (const Matrix< T > &D, const std::vector< int > &N, const std::vector< T > &Z, McRng &rng)
	Reference defaults: 1e5 samples, the IMCI proposal.

Detailed Description

Monte Carlo Integration estimate of the normalizing constant of a closed product-form network (Ross, Wang and Yao; MonteQueue 2.0).

Templated port of matlab/src/api/pfqn/pfqn_mci.m, cross-checked against jar/src/main/java/jline/api/pfqn/nc/Pfqn_mci.java.

The McKenna-Mitra integral form of the constant,

G(N) = 1/prod_r N_r! int_{R_+^M} e^{-sum_i v_i} prod_r (sum_i v_i D(i,r) + Z_r)^{N_r} dv,

is estimated by importance sampling with an independent exponential proposal of rate gamma_i per station. Writing V ~ prod_i Exp(gamma_i), one draw contributes

lZ = -sum_i (1 - gamma_i) V_i - sum_i log gamma_i - sum_r log N_r!

• sum_r N_r log(sum_i V_i D(i,r) + Z_r),

and lG = log mean exp(lZ). The proposal rates come from a Bard-Schweitzer pre-solve: gamma_i = max(0.01, 1 - U_i) in the IMCI variant, and the saturation-aware gamma_i = 1/sqrt(max N) when U_i > 0.9 in the plain MCI variant. Overshooting rates make the estimator heavy-tailed, which is why the reference clamps them.

Deviations from the reference, both deliberate:

• MATLAB caches the uniform matrix in a PERSISTENT variable VL and reuses it across calls (slicing VL(1:I,1:M)), so two successive calls with the same shape return the SAME estimate and a call with a larger I silently reuses the old columns. That is global hidden state; here the deviates come from the caller's generator, one draw per station per sample.
• The 'rm' variant is REFUSED for M > 1 rather than reproduced. Its rate line, tput = N./(sum(D,1)+Z+max(D,1)*(sum(N)-1)), broadcasts a (1 x R) row against the (M x R) matrix max(D,1) and yields an (M x R) "throughput" whose product D*tput' is (M x M); the subsequent loop then reads gamma off the wrong axis. The expression is only dimensionally meaningful for a single station, which is the repairman model the variant is named for, so that is the only shape accepted here. See the report accompanying this port; nothing has been substituted for the M > 1 case.

Arithmetic: INEXACT BY CONSTRUCTION. The value is a random variable, the integrand is evaluated through logs, and the proposal rates come from a tolerance-stopped Bard-Schweitzer solve; all three require transcendental arithmetic.

RNG contract: see pfqn_mc_common.h. Comparable to MATLAB only in distribution, never stream for stream; reproducible within this port only when the generator is passed in the same state.

Definition in file pfqn_mci.h.