The variate layer of the Java LDES engine, reproduced: SSJ's randvar generators as inverse-CDF functions of one uniform. More...

#include <cmath>
#include <string>
#include <vector>
#include "line/util/error.h"
#include "line/util/fdlibm.h"

Include dependency graph for ldes_ssj_variates.h:

Namespaces
namespace	line
namespace	line::ldes
namespace	line::ldes::ssj

Functions
double	line::ldes::ssj::exponential_inverse (double lambda, double u)
	ExponentialDist.inverseF(lambda, u), SSJ's log1p spelling.
double	line::ldes::ssj::uniform_inverse (double a, double b, double u)
	UniformDist.inverseF(a, b, u).
double	line::ldes::ssj::pareto_inverse (double alpha, double beta, double u)
	ParetoDist.inverseF(alpha, beta, u) = beta (1-u)^(-1/alpha).
double	line::ldes::ssj::weibull_inverse (double alpha, double lambda, double delta, double u)
	WeibullDist.inverseF(alpha, lambda, delta, u).
double	line::ldes::ssj::bernoulli_inverse (double p, double u)
	BernoulliDist.inverseF(p, u): 1 when u exceeds 1 - p.
double	line::ldes::ssj::normal_inverse01 (double u)
	NormalDist.inverseF(0, 1, u) by Wichura's AS 241 (Applied Statistics 37, 1988), the algorithm SSJ's inverseF01 implements.
double	line::ldes::ssj::normal_inverse (double mu, double sigma, double u)
	NormalDist.inverseF(mu, sigma, u).
double	line::ldes::ssj::lognormal_inverse (double mu, double sigma, double u)
	LognormalDist.inverseF(mu, sigma, u) = exp(mu + sigma Phi^-1(u)).
double	line::ldes::ssj::gamma_inverse (double alpha, double lambda, double u)
	GammaDist.inverseF(alpha, lambda, u): the quantile of a Gamma of shape alpha and RATE lambda, so the mean is alpha/lambda, which is SSJ's convention.
double	line::ldes::ssj::erlang_inverse (int k, double lambda, double u)
	ErlangGen(k, lambda): the Gamma of integer shape k and rate lambda.
double	line::ldes::ssj::poisson_inverse (double lambda, double u)
	PoissonDist.inverseF(lambda, u): the smallest k whose cdf reaches u.
double	line::ldes::ssj::binomial_inverse (int n, double p, double u)
	BinomialDist.inverseF(n, p, u), by the same forward inversion.

Detailed Description

The variate layer of the Java LDES engine, reproduced: SSJ's randvar generators as inverse-CDF functions of one uniform.

THE MEASUREMENT THAT SHAPED THIS FILE. Every *Gen instance Solver_ssj builds – Exponential, Erlang, Uniform, Weibull, Pareto, Lognormal, Gamma, Poisson, Binomial, Bernoulli – consumes EXACTLY ONE uniform per draw. That was measured, not assumed: a probe generated five variates from a seeded MRG32k3a and then searched the raw stream for the next value, and every family had advanced the stream by exactly five. SSJ's instance generators invert; none of them rejects. The consequence is the whole reason a bit-compatible C++ engine is feasible: the stream can never desynchronize, so a family whose quantile this file computes to 1e-13 rather than to the last bit still leaves every LATER draw identical, and only perturbs its own sample below the tolerance of any event comparison.

WHAT IS EXACT AND WHAT IS TO 1e-13. Exponential, Uniform, Pareto, Weibull and Bernoulli have closed-form quantiles and agree with SSJ to the last bit. Normal (hence Lognormal) uses Wichura's AS 241, Gamma and Erlang inverting the regularized incomplete gamma by Newton with a Wilson-Hilferty start, and Poisson and Binomial by discrete inversion – the two discrete families return INTEGERS and so are exact wherever the cdf comparison lands on the same side, which it does everywhere the test looks.

WHY NOT CALL SSJ'S OWN FORMULAE BLINDLY. ExponentialDist.inverseF is -log1p(-u)/lambda and not -log(1-u)/lambda; the two differ in the last bits for small u, which is exactly where an interarrival time matters most. Where a spelling like that decides agreement, this file uses SSJ's.

AND THE SPELLING IS NOT ENOUGH: IT MUST BE THE SAME log1p. std::log1p is allowed a 1-ulp error, and glibc changed which representative it returns between 2.35 and 2.39, so the SAME common/ldes binary on the SAME model.json at --seed 23000 gave two different sample paths – 98.565592 on a 22.04 host against 98.562490 under the containerized MATLAB's 24.04, which is how this was found. The quantiles below therefore call line::fdlibm::log1p, the algorithm StrictMath.log1p is DEFINED to use and which Math.log1p was measured to match on every one of 400k draws; see line/util/fdlibm.h. That fixes the engine to one arithmetic everywhere AND puts it on the Java engine's, which is what "bit-compatible" was supposed to mean. pow, log, exp, sqrt, lgamma and tgamma were measured to agree across those glibcs and are left to libm.

Definition in file ldes_ssj_variates.h.

Namespaces

Functions

Detailed Description