Class TruncatedGaussian

Class Documentation

class TruncatedGaussian

Represents a multidimensional Gaussian function truncated at zero.

This is typically used to represent the posterior probability of amplitude parameters, given a flat prior; we require that the amplitudes each be positive but otherwise do not modify the Gaussian likelihood.

Currently only 1 and 2 dimensions are supported, and all dimensions must be truncated. Computing integrals is the only operation for which > 2-d is not implemented, but the integrals must be computed upon construction, so we can’t support any other operations for > 2-d either.

Many operations on TruncatedGaussians are defined in -log space, as underflow/overflow problems will often occur in the non-log forms.

See modelfitTruncatedGaussianMath for implementation notes

Public Types

enum SampleStrategy

Enum that describes different ways of sampling from a multidimensional TruncatedGaussian.

Values:

DIRECT_WITH_REJECTION

Draw from the untruncated Gaussian, and discard negative draws. This approach is fast for minimally-truncated Gaussians, and produces ideally distributed samples with unit weights (because we discard with replacement). This approach can be extremely inefficient when the untruncated fraction is small, and cannot be used when the Hessian matrix is singular (as this sets the untruncated fraction to zero).

ALIGN_AND_WEIGHT

Create a similar Gaussian with no x-y covariance, and importance sample by drawing from the independent 1-d truncated Gaussians. This approach is faster for Gaussians that are heavily truncated, but it requires the use of per-sample weights and may produce “noisier” samples (in that regions with low weight may have many samples, while regions with high weight have few). It may also be slower for minimally truncated Gaussians because each 1-d variate requires an inverse erfc evaluation. In 1-d, this produces unit weights and ideally distributed samples, but still uses an inverse erfc approach rather than rejection to create samples.

typedef TruncatedGaussianSampler Sampler
typedef TruncatedGaussianLogEvaluator LogEvaluator
typedef TruncatedGaussianEvaluator Evaluator

Public Functions

TruncatedGaussian::Sampler sample(SampleStrategy strategy) const

Create a Sampler object that uses the given strategy.

Warning

do not call this method with DIRECT_WITH_REJECTION unless you are certain the Gaussian is far from degenerate, as an endless loop may result in this case. It is generally safer to use the overload that selects the strategy automatically.

TruncatedGaussian::Sampler sample(Scalar minRejectionEfficiency = 0.1) const

Create a Sampler object that determines the strategy to use automatically.

If the efficiency (number of accepted samples divided by number of total samples) is would (on average) be greater than the given value, DIRECT_WITH_REJECTION is used. If not, ALIGN_AND_WEIGHT is used. Note that the sampling effeciency for DIRECT_WITH_REJECTION is exactly the same as the “untruncated fraction”.

TruncatedGaussian::LogEvaluator evaluateLog() const

Create a LogEvaluator object that can be used to efficiently evaluate the -log of the function.

TruncatedGaussian::Evaluator evaluate() const

Create an Evaluator object that can be used to efficiently evaluate the function.

int getDim() const

Return the dimensionality of the function.

Vector maximize() const

Return the location of the maximum of the truncated Gaussian.

This is simply the untruncated location parameter mu if all of its elements are positive; otherwise, one or more elements will be zero (and the rest may not be same as the elements of mu).

Scalar getUntruncatedFraction() const

Return the fraction of the Gaussian integral that was truncated by the bounds.

In series form (see fromSeriesParameters), this is

\[ \frac{\int_0^{\infty} e^{-q(\alpha)} d\alpha}{\int_{-infty}^{\infty} e^{-q(\alpha)} d\alpha} \]
while in standard parameter form it is simply the normalization \(k\) described in fromStandardParameters().

Scalar getLogPeakAmplitude() const

Return the -log of the peak amplitude of the untruncated function.

In series form, this is equal to \(q(\mu)\) where \(\mu=-H^{-1}g\). In standard parameter form, this is \(\frac{1}{2}\ln\left|2\pi\Sigma\right| - \ln k\).

Scalar getLogIntegral() const

Return the -log of the integral of the truncated function.

More precisely, this is simply

\[ -\ln\int_{-\infty}^{\infty} f(\alpha) d\alpha \]
In series form, this is equivalent to
\[ -\ln\int_0^{\infty} e^{-q(\alpha)} d\alpha \]
In standard parameter form it is 0 by construction.

~TruncatedGaussian()

Public Static Functions

static TruncatedGaussian fromSeriesParameters(Scalar q0, Vector const &gradient, Matrix const &hessian)

Create from the first and second logarithmic derivatives of the Gaussian.

This causes the TruncatedGaussian to model the function:

\[ f(\alpha) = 1_{\alpha \ge 0}(\alpha) \; e^{-q(0) -g^T \alpha -\frac{1}{2}\alpha^T H \alpha} = 1_{\alpha \ge 0}(\alpha) \; e^{-q(\alpha)} \]
Where \(1_{\alpha \ge 0}(\alpha)\) is the indicator function
\[\begin{split} 1_{\alpha \ge 0}(\alpha) \equiv \begin{cases} 1 & \text{if $\alpha \ge 0$}\\ 0 & \text{if $\alpha \lt 0$} \end{cases} \end{split}\]
Note that this implies
\[ \left.\frac{\partial q}{\partial \alpha} \right|_{\alpha=0} = g \]
\[ \left.\frac{\partial^2 q}{\partial \alpha^2} \right|_{\alpha=0} = H \]
with (for \(\alpha \ge 0\))
\[ q(\alpha) \equiv -\ln f(\alpha) \]

Note that unlike that constructed by fromStandardParameters(), in this case \(f(\alpha)\) is NOT normalized to integrate to one.

Parameters

  • [in] q0: Zeroth-order term in the expansion

  • [in] gradient: Vector of first derivatives \(g\)

  • [in] hessian: Matrix of second derivatives \(H\) (symmetric)

static TruncatedGaussian fromStandardParameters(Vector const &mean, Matrix const &covariance)

Create from the “standard” mean and covariance parameters of the normal distribution.

This causes the TruncatedGaussian to model the function:

\[ f(\alpha) = \frac{1_{\alpha \ge 0}(\alpha)}{k\left|2\pi\Sigma\right|^{1/2}} e^{-\frac{1}{2}(\alpha-\mu)^T \Sigma^{-1} (\alpha-\mu)} = \frac{1_{\alpha \ge 0}(\alpha)}{k}\;\Phi(\alpha-\mu,\Sigma) \]
Where \(1_{\alpha \ge 0}(\alpha)\) is the indicator function
\[\begin{split} 1_{\alpha \ge 0}(\alpha) \equiv \begin{cases} 1 & \text{if $\alpha \ge 0$}\\ 0 & \text{if $\alpha \lt 0$} \end{cases} \end{split}\]
and \(k\) is the normalization constant that \(f(\alpha)\) integrates to unity:
\[ k \equiv \int_0^{\infty}\!\!\Phi(\alpha-\mu,\Sigma) d\alpha \]

Note that unlike that constructed by fromSeriesParameters(), in this case \(f(\alpha\) is normalized to integrate to one.

Parameters

  • [in] mean: Mean vector \(\mu\)

  • [in] covariance: Covariance matrix \(\Sigma\) (symmetric)