File TruncatedGaussian.h¶
-
namespace
lsst Class for a simple mapping implementing a generic AstrometryTransform.
Remove all non-astronomical counts from the Chunk Exposure’s pixels.
Forward declarations for lsst::utils::Cache
For details on the Cache class, see the Cache.h file.
It uses a template rather than a pointer so that the derived classes can use the specifics of the transform. The class simplePolyMapping overloads a few routines.
A base class for image defects
Numeric constants used by the Integrate.h integrator routines.
Compute Image Statistics
- Note
Gauss-Kronrod-Patterson quadrature coefficients for use in quadpack routine qng. These coefficients were calculated with 101 decimal digit arithmetic by L. W. Fullerton, Bell Labs, Nov 1981.
- Note
The Statistics class itself can only handle lsst::afw::image::MaskedImage() types. The philosophy has been to handle other types by making them look like lsst::afw::image::MaskedImage() and reusing that code. Users should have no need to instantiate a Statistics object directly, but should use the overloaded makeStatistics() factory functions.
-
namespace
meas -
namespace
modelfit -
class
TruncatedGaussian - #include <TruncatedGaussian.h>
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)
Private Functions
-
lsst::meas::modelfit::TruncatedGaussian::TruncatedGaussian(PTR(Impl) impl)
-
PTR(Impl)¶
Friends
-
friend
lsst::meas::modelfit::TruncatedGaussianSampler
-
friend
lsst::meas::modelfit::TruncatedGaussianLogEvaluator
-
enum
-
class
TruncatedGaussianEvaluator - #include <TruncatedGaussian.h>
Helper class for evaluating the -log of a TruncatedGaussian.
Public Functions
-
TruncatedGaussianEvaluator(TruncatedGaussian const &parent)
Private Members
-
TruncatedGaussianLogEvaluator
_internal¶
-
-
class
TruncatedGaussianLogEvaluator - #include <TruncatedGaussian.h>
Helper class for evaluating the -log of a TruncatedGaussian.
Public Functions
-
TruncatedGaussianLogEvaluator(TruncatedGaussian const &parent)
-
-
class
TruncatedGaussianSampler - #include <TruncatedGaussian.h>
Helper class for drawing samples from a TruncatedGaussian.
Public Functions
-
TruncatedGaussianSampler(TruncatedGaussian const &parent, TruncatedGaussian::SampleStrategy strategy)
-
Scalar
operator()(afw::math::Random &rng, ndarray::Array<Scalar, 1, 1> const &alpha) const Draw a single sample from a TruncatedGaussian.
- Return
the weight of the sample (always betweeen 0 and 1)
- Parameters
[in] rng: Random number generator[out] alpha: Output sample vector to fill
-
void
operator()(afw::math::Random &rng, ndarray::Array<Scalar, 2, 1> const &alpha, ndarray::Array<Scalar, 1, 1> const &weights, bool multiplyWeights = false) const Draw multiple samples from a TruncatedGaussian.
- Parameters
[in] rng: Random number generator[out] alpha: Output sample vector to fill; first dimension sets the number of samples[out] weights: Output weight vector to fill[in] multiplyWeights: If true, multiply the weights vector by the weights rather than fill it.
-
~TruncatedGaussianSampler()
Private Functions
-
PTR(Impl)¶
-
-
class
-
namespace