Chen-O'Cinneide REGULARIZATION: a Markov chain Monte Carlo estimator of the class throughputs X(r) = G(N-e_r)/G(N) and of the mean queue lengths Q(i,r) of a CLOSED multiclass product-form (BCMP, no type changes) network. More...

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

Classes
struct	line::pfqn::McmcResult< T >
	Estimates of pfqn_mcmc together with their batch-means intervals. More...

Namespaces
namespace	line
namespace	line::pfqn

Functions
template<class T>
McmcResult< T >	line::pfqn::pfqn_mcmc (const Matrix< T > &L, const std::vector< int > &N, const std::vector< T > &Z, const std::vector< double > &s, std::size_t samples, std::size_t nbatches, double burnin, McRng &rng)
	Chen-O'Cinneide REGULARIZATION: a Markov chain Monte Carlo estimator of the class throughputs X(r) = G(N-e_r)/G(N) and of the mean queue lengths Q(i,r) of a CLOSED multiclass product-form (BCMP, no type changes) network.

Variables
constexpr std::size_t	line::pfqn::MCMC_DEFAULT_BATCHES = 30
	Schmeiser (1982), the batch count used in the tables of the paper.
constexpr double	line::pfqn::MCMC_DEFAULT_BURNIN = 0.1
	Warm-up fraction discarded before accumulation starts.

Detailed Description

Chen-O'Cinneide REGULARIZATION: a Markov chain Monte Carlo estimator of the class throughputs X(r) = G(N-e_r)/G(N) and of the mean queue lengths Q(i,r) of a CLOSED multiclass product-form (BCMP, no type changes) network.

W. Chen, C. A. O'Cinneide, "Towards a Polynomial-Time Randomized Algorithm for Closed Product-Form Networks", ACM TOMACS 8(3):227-253, 1998.

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

The three steps of the paper are:

I. CONSTRUCT THE REGULARIZED NETWORK. Write rho(i,r) for the surrogate traffic intensity of class r at station i – here the service demand, since rho = lambda/mu is a visit ratio over a service rate – and rho(r) = sum_i rho(i,r). The regularized network has the same stations, classes and populations, UNIT service rates at every station, the processor-sharing discipline, and a routing matrix that depends on the destination only,

P*(i->m | class r) = rho(m,r)/rho(r).

By Theorem 2.1 it is a REVERSIBLE chain with the SAME steady-state distribution as the original network, and its throughputs satisfy Theta*(r) = rho(r)*Theta(r).

II. SIMULATE IT at service-completion epochs. With Y(i,r) the number of class-r jobs at station i, Y(i) their total and Psi_i(k) = min(s_i,k) the number of busy servers,

r(i,r) = Y(i,r)/Y(i) * Psi_i(Y(i)),   r(r) = sum_i r(i,r),
r      = sum_i Psi_i(Y(i)),

the next completion is of class r at station i with probability r(i,r)/r, and the conditional expected time to it is 1/r. Equation (10) of the paper is the holding-time weighted ratio estimator

Theta*(r) = sum_t r(r,t)/r(t)  /  sum_t 1/r(t),

and the same weights give the time-average queue lengths, which need no transformation at all because the two networks share their steady state.

III. TRANSFORM BACK: X(r) = Theta*(r)/rho(r).

Because P* forgets the station of origin and every station serves at unit rate, the regularized chain has neither the slowly mixing routing chain nor the customer-trapping slow station that make the original chain converge slowly. The paper proves O(N^2*M^3) mixing in two special cases (Section 4) and reports the general behaviour experimentally (Section 5).

Delay (infinite-server) demand enters as ONE extra station with s = infinity and demand Z. Aggregating infinite-server stations that way is exact in the product form, since their joint term is multinomial in the per-class totals.

THERE IS NO NORMALIZING CONSTANT HERE, and that is the single most important thing to know about the method: the estimator is a ratio, G itself never appears. pfqn_nc's Mcmc arm supplies an BLE lG alongside, which is not part of the paper and cancels out of every mean value.

Confidence: the run is split into non-overlapping batches (Schmeiser 1982, 30 by default, the count used in the tables of the paper), the batch means of the ratio estimator give a standard error, and the intervals are the paper's two-sigma ones. The estimator is a ratio of correlated averages, so it carries an O(1/samples) bias on top of the initialization bias; the paper ignores both, this port additionally discards a warm-up fraction (10% by default).

Arithmetic: INEXACT BY CONSTRUCTION. The value is a random variable and the batch-means interval needs a square root, so the estimator is meaningless – not merely inaccurate – in an exact field, exactly as the rest of the pfqn_nc estimator ladder is.

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

Classes

Namespaces

Functions

Variables

Detailed Description