Sampled McKenna-Mitra integral form of the normalizing constant of a repairman (single-queue plus delay) model. More...

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

Namespaces
namespace	line
namespace	line::pfqn

Functions
template<class T>
NcResult< T >	line::pfqn::pfqn_mmsample2 (const Matrix< T > &L, const std::vector< int > &N, const std::vector< T > &Z, std::size_t samples, McRng &rng)
	Sampled McKenna-Mitra integral form of the normalizing constant of a repairman (single-queue plus delay) model.
template<class T>
NcResult< T >	line::pfqn::pfqn_mmsample2 (const Matrix< T > &L, const std::vector< int > &N, const std::vector< T > &Z, McRng &rng)
	Reference call shape with an explicit grid size.

Detailed Description

Sampled McKenna-Mitra integral form of the normalizing constant of a repairman (single-queue plus delay) model.

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

For a single queue the McKenna-Mitra form collapses to the one-dimensional integral

G(N) = 1/prod_r N_r! int_0^inf e^{-v} prod_r ((Z_r + L_r) v)^{N_r} dv,

whose log-integrand at the grid point v is

f(v) = -v + sum(N) log v + sum_r N_r log(Z_r + L_r).

The reference lays a grid of ceil(samples/2) uniform points on [0,1) plus ceil(samples/2) points log-spaced over [1, 1e5], rescaling the demands so that every coefficient is at least one, and reports

lG = max_i (du_i + f(v_i)) - sum_r log N_r! + sum(N) log(scaleFactor).

DEFECT IN THE REFERENCE, reproduced here rather than repaired. The estimator takes the MAXIMUM of the weighted log-integrand instead of summing the quadrature contributions log sum_i du_i exp(f(v_i)). Two consequences:

the width of the integrand is discarded, so the estimate is short by the Laplace factor sqrt(2 pi sum(N)): with the exact value being log((sum N)!) + sum_r N_r log(Z_r+L_r) - sum_r log N_r!, the returned value is Stirling's leading term without the sqrt(2 pi n) correction, a systematic relative shortfall of about 1/sqrt(2 pi sum(N)) in the log and a factor sqrt(2 pi sum(N)) in G itself;
the grid v is NOT sorted (uniform points first, then an increasing log-spaced block), so du = [0, diff(v)] alternates in sign and is not a quadrature weight at all; it merely perturbs the argmax.

Because of this the routine does not converge to the exact constant as the sample count grows, and its test in this tree asserts the reference's own behaviour and the size of the shortfall, not agreement with pfqn_ca.

Only the FIRST row of L is read, matching the reference: this is a single-queue model and additional rows are ignored.

Arithmetic: INEXACT BY CONSTRUCTION. The grid is random, the integrand is evaluated in the log domain, and the log-spaced block is generated by exponentiation.

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_mmsample2.h.

Namespaces

Functions

Detailed Description