fluid_mvn_rectangle.h File Reference

Port of fluid_mvn_rectangle.m: the rectangle probability P(a <= Y <= b) for Y ~ Normal(m, C), the cell integral behind getProbAggr under the moment-closure methods. More...

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <vector>
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for fluid_mvn_rectangle.h:

Go to the source code of this file.

Namespaces
namespace	line
namespace	line::fluid

Functions
double	line::fluid::fluid_mvn_rectangle (const std::vector< double > &m, const Matrix< double > &C, const std::vector< double > &a, const std::vector< double > &b, std::size_t npoints=FLUID_MVN_POINTS)
	P(a <= Y <= b) for Y ~ Normal(m, C).

Variables
const std::size_t	line::fluid::FLUID_MVN_POINTS = 4096
	Lattice points per antithetic pair.

Detailed Description

Port of fluid_mvn_rectangle.m: the rectangle probability P(a <= Y <= b) for Y ~ Normal(m, C), the cell integral behind getProbAggr under the moment-closure methods.

The integral has no closed form beyond one dimension, so it is evaluated by the separation-of-variables transformation of Genz (1992): the Cholesky factor of C turns the rectangle into an iterated integral over the unit cube whose integrand is a product of normal-CDF differences, and the first coordinate is integrated exactly. The remaining cube is integrated with a DETERMINISTIC Richtmyer lattice rule, frac(k sqrt(p_j)) over the first primes, averaged with its antithetic reflection. Determinism is required here, not merely convenient: the four codebases must return the same number, and a randomized rule would make them agree only in distribution.

C may be SINGULAR, which is the common case: a closed population fixes the sum of the station coordinates, so the covariance of a station holding a whole class is rank deficient. A coordinate whose CONDITIONAL variance vanishes is not integrated; it is a hard constraint, contributing 1 when the conditional mean falls inside its interval and 0 otherwise.

Definition in file fluid_mvn_rectangle.h.