A class for handling hybrid (transformed) Gaussians with bounds. More...

#include <InvLogitGaussianJointPdf.h>

Inheritance diagram for QUESO::InvLogitGaussianJointPdf< V, M >:

Collaboration diagram for QUESO::InvLogitGaussianJointPdf< V, M >:

Public Member Functions
virtual void	print (std::ostream &os) const
	Prints the distribution. More...

Constructor/Destructor methods
	InvLogitGaussianJointPdf (const char *prefix, const BoxSubset< V, M > &domainBoxSubset, const V &lawExpVector, const V &lawVarVector)
	Constructor. More...

	InvLogitGaussianJointPdf (const char *prefix, const BoxSubset< V, M > &domainBoxSubset, const V &lawExpVector, const M &lawCovMatrix)
	Constructor. More...

	~InvLogitGaussianJointPdf ()
	Destructor. More...

Math methods
double	actualValue (const V &domainVector, const V domainDirection, V gradVector, M hessianMatrix, V hessianEffect) const
	Actual value of the (transformed) Gaussian PDF. More...

double	lnValue (const V &domainVector, const V domainDirection, V gradVector, M hessianMatrix, V hessianEffect) const
	Logarithm of the value of the (transformed) Gaussian PDF (scalar function). More...

virtual void	distributionMean (V &meanVector) const
	Mean value of the underlying random variable. More...

virtual void	distributionVariance (M &covMatrix) const
	Covariance matrix of the underlying random variable. More...

double	computeLogOfNormalizationFactor (unsigned int numSamples, bool updateFactorInternally) const
	Computes the logarithm of the normalization factor. More...

void	updateLawExpVector (const V &newLawExpVector)
	Updates the mean of the Gaussian (not transformed) with the new value `newLawExpVector`. More...

void	updateLawCovMatrix (const M &newLawCovMatrix)
	Updates the lower triangular matrix from Cholesky decomposition of the covariance matrix to the new value `newLowerCholLawCovMatrix`. More...

const M &	lawCovMatrix () const
	Returns the covariance matrix; access to protected attribute m_lawCovMatrix. More...

const V &	lawExpVector () const
	Access to the vector of mean values of the Gaussian (not transformed) and private attribute: m_lawExpVector. More...

const V &	lawVarVector () const
	Access to the vector of variance values and private attribute: m_lawVarVector. More...

Public Member Functions inherited from QUESO::BaseJointPdf< V, M >
	BaseJointPdf (const char *prefix, const VectorSet< V, M > &domainSet)
	Default constructor. More...

virtual	~BaseJointPdf ()
	Destructor. More...

virtual void	setNormalizationStyle (unsigned int value) const
	Sets a value to be used in the normalization style (stored in the protected attribute m_normalizationStyle.) More...

void	setLogOfNormalizationFactor (double value) const
	Sets a logarithmic value to be used in the normalization factor (stored in the protected attribute m_normalizationStyle.) More...

Public Member Functions inherited from QUESO::BaseScalarFunction< V, M >
	BaseScalarFunction (const char *prefix, const VectorSet< V, M > &domainSet)
	Default constructor. More...

virtual	~BaseScalarFunction ()
	Destructor. More...

const VectorSet< V, M > &	domainSet () const
	Access to the protected attribute `m_domainSet:` domain set of the scalar function. More...

Private Attributes
V *	m_lawExpVector

V *	m_lawVarVector

bool	m_diagonalCovMatrix

const M *	m_lawCovMatrix

const BoxSubset< V, M > &	m_domainBoxSubset

Additional Inherited Members
Protected Member Functions inherited from QUESO::BaseJointPdf< V, M >
double	commonComputeLogOfNormalizationFactor (unsigned int numSamples, bool updateFactorInternally) const
	Common method (to the derived classes) to compute the logarithm of the normalization factor. More...

Protected Attributes inherited from QUESO::BaseJointPdf< V, M >
unsigned int	m_normalizationStyle

double	m_logOfNormalizationFactor

Protected Attributes inherited from QUESO::BaseScalarFunction< V, M >
const BaseEnvironment &	m_env

std::string	m_prefix

const VectorSet< V, M > &	m_domainSet
	Domain set of the scalar function. More...

Detailed Description

template<class V = GslVector, class M = GslMatrix>
class QUESO::InvLogitGaussianJointPdf< V, M >

A class for handling hybrid (transformed) Gaussians with bounds.

This class allows the mathematical definition of a Gaussian PDF in dimensions of a VectorSubset that are unbounded, and a logit transformed Gaussian PDF in dimensions that have a lower or upper (or both) bound. For the bounded directions, the PDF is

$p(x) = \frac{1}{2 \sigma^2 \pi} \exp \left( - \frac12 (\mbox{logit} x - \mu)^2 \right) \mbox{logit}'(x), \quad x \in [a, b],$

where logit is defined as

$\mbox{logit}(x) = \log \left( \frac{x - a}{b - x} \right), \quad x \in [a, b].$

Definition at line 58 of file InvLogitGaussianJointPdf.h.

Constructor & Destructor Documentation

template<class V , class M >

QUESO::InvLogitGaussianJointPdf< V, M >::InvLogitGaussianJointPdf	(	const char *	prefix,
		const BoxSubset< V, M > &	domainBoxSubset,
		const V &	lawExpVector,
		const V &	lawVarVector
	)

Constructor.

Constructs a new object, given a prefix and the domain of the PDF, a vector of mean values, lawExpVector (for the Gaussian, not the transformed Gaussian), and a vector of covariance values lawVarVector (an alternative representation for a diagonal covariance matrix).

Definition at line 33 of file InvLogitGaussianJointPdf.C.

   :
   BaseJointPdf<V,M>(((std::string)(prefix)+"invlogit_gau").c_str(),
       domainBoxSubset),
   m_lawExpVector(new V(lawExpVector)),
   m_lawVarVector(new V(lawVarVector)),
   m_diagonalCovMatrix(true),
   m_lawCovMatrix(m_domainSet.vectorSpace().newDiagMatrix(lawVarVector)),
   m_domainBoxSubset(domainBoxSubset)
 {
 }

template<class V , class M >

QUESO::InvLogitGaussianJointPdf< V, M >::InvLogitGaussianJointPdf	(	const char *	prefix,
		const BoxSubset< V, M > &	domainBoxSubset,
		const V &	lawExpVector,
		const M &	lawCovMatrix
	)

Constructor.

Constructs a new object, given a prefix and the image set of the vector realizer, a vector of mean values, lawExpVector (for the Gaussian, not the transformed Gaussian), and a covariance matrix, lawCovMatrix.

Definition at line 51 of file InvLogitGaussianJointPdf.C.

   :
   BaseJointPdf<V,M>(((std::string)(prefix)+"invlogit_gau").c_str(),
       domainBoxSubset),
   m_lawExpVector(new V(lawExpVector)),
   m_lawVarVector(domainBoxSubset.vectorSpace().newVector(INFINITY)), // FIX ME
   m_diagonalCovMatrix(false),
   m_lawCovMatrix(new M(lawCovMatrix)),
   m_domainBoxSubset(domainBoxSubset)
 {
 }

template<class V , class M >

QUESO::InvLogitGaussianJointPdf< V, M >::~InvLogitGaussianJointPdf ( )

Destructor.

Definition at line 68 of file InvLogitGaussianJointPdf.C.

 {
   delete m_lawCovMatrix;
   delete m_lawVarVector;
   delete m_lawExpVector;
 }

Member Function Documentation

template<class V , class M >

double QUESO::InvLogitGaussianJointPdf< V, M >::actualValue	(	const V &	domainVector,
		const V *	domainDirection,
		V *	gradVector,
		M *	hessianMatrix,
		V *	hessianEffect
	)		const

virtual

Actual value of the (transformed) Gaussian PDF.

This method calls lnValue() and applies the exponential to it.

Implements QUESO::BaseJointPdf< V, M >.

Definition at line 117 of file InvLogitGaussianJointPdf.C.

 {
   double returnValue;
 
   returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
 
   return returnValue;
 }

template<class V , class M >

double QUESO::InvLogitGaussianJointPdf< V, M >::computeLogOfNormalizationFactor	(	unsigned int	numSamples,
		bool	updateFactorInternally
	)		const

virtual

Computes the logarithm of the normalization factor.

This routine calls BaseJointPdf::commonComputeLogOfNormalizationFactor().

Implements QUESO::BaseJointPdf< V, M >.

Definition at line 257 of file InvLogitGaussianJointPdf.C.

References QUESO::BaseJointPdf< V, M >::commonComputeLogOfNormalizationFactor().

 {
   double value = 0.;
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
     *m_env.subDisplayFile() << "Entering GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
                             << std::endl;
   }
   value = BaseJointPdf<V,M>::commonComputeLogOfNormalizationFactor(numSamples, updateFactorInternally);
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
     *m_env.subDisplayFile() << "Leaving GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
                             << ", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
                             << std::endl;
   }
 
   return value;
 }

template<class V , class M >

void QUESO::InvLogitGaussianJointPdf< V, M >::distributionMean ( V & meanVector ) const

virtual

Mean value of the underlying random variable.

Reimplemented from QUESO::BaseJointPdf< V, M >.

Definition at line 233 of file InvLogitGaussianJointPdf.C.

References queso_not_implemented.

 {
   // AFAIK there's no simple closed form here, and taking the inverse
   // transformation of the mean in the transformed space probably
   // isn't accurate enough in cases where the mean is too near the
   // bounds.
   queso_not_implemented();
 }

template<class V , class M >

void QUESO::InvLogitGaussianJointPdf< V, M >::distributionVariance ( M & covMatrix ) const

virtual

Covariance matrix of the underlying random variable.

Reimplemented from QUESO::BaseJointPdf< V, M >.

Definition at line 245 of file InvLogitGaussianJointPdf.C.

References queso_not_implemented.

 {
   // AFAIK there's no simple closed form here, and taking the inverse
   // transformation of the variance in the transformed space probably
   // isn't accurate enough in cases where the mean is too near the
   // bounds.
   queso_not_implemented();
 }

template<class V , class M >

const M & QUESO::InvLogitGaussianJointPdf< V, M >::lawCovMatrix ( ) const

Returns the covariance matrix; access to protected attribute m_lawCovMatrix.

Definition at line 295 of file InvLogitGaussianJointPdf.C.

Referenced by QUESO::TransformedScaledCovMatrixTKGroup< V, M >::setPreComputingPosition().

 {
   return *m_lawCovMatrix;
 }

template<class V , class M >

const V & QUESO::InvLogitGaussianJointPdf< V, M >::lawExpVector ( ) const

Access to the vector of mean values of the Gaussian (not transformed) and private attribute: m_lawExpVector.

Definition at line 77 of file InvLogitGaussianJointPdf.C.

 {
   return *m_lawExpVector;
 }

template<class V , class M >

const V & QUESO::InvLogitGaussianJointPdf< V, M >::lawVarVector ( ) const

Access to the vector of variance values and private attribute: m_lawVarVector.

Definition at line 84 of file InvLogitGaussianJointPdf.C.

 {
   return *m_lawVarVector;
 }

template<class V , class M >

double QUESO::InvLogitGaussianJointPdf< V, M >::lnValue	(	const V &	domainVector,
		const V *	domainDirection,
		V *	gradVector,
		M *	hessianMatrix,
		V *	hessianEffect
	)		const

virtual

Logarithm of the value of the (transformed) Gaussian PDF (scalar function).

The ln(value) comes from a summation of the Gaussian density:

$lnValue =- \sum_i \frac{1}{\sqrt{|covMatrix|} \sqrt{2 \pi}} exp(-\frac{(domainVector_i - lawExpVector_i) * covMatrix^{-1} * (domainVector_i - lawExpVector_i) }{2} + \log \mbox{logit}'(domainVector_i),$

where the $ covMatrix $ may recovered via this->lawVarVector(), in case of diagonal matrices or via this->m_lawCovMatrix, otherwise.

Implements QUESO::BaseJointPdf< V, M >.

Definition at line 133 of file InvLogitGaussianJointPdf.C.

References QUESO::queso_isfinite().

 {
   double returnValue;
   double lnDeterminant = 0.0;
   V transformedDomainVector(domainVector);
 
   V min_domain_bounds(this->m_domainBoxSubset.minValues());
   V max_domain_bounds(this->m_domainBoxSubset.maxValues());
 
   double lnjacobian = 0.0;
   for (unsigned int i = 0; i < domainVector.sizeLocal(); i++) {
     double min_val = min_domain_bounds[i];
     double max_val = max_domain_bounds[i];
 
     if (queso_isfinite(min_val) &&
         queso_isfinite(max_val)) {
 
       if (domainVector[i] == min_val || domainVector[i] == max_val) {
         // Exit early if we can
         return -INFINITY;
       }
 
         // Left- and right-hand sides are finite.  Do full transform.
         transformedDomainVector[i] = std::log(domainVector[i] - min_val) -
             std::log(max_val - domainVector[i]);
 
         lnjacobian += std::log(max_val - min_val) -
           std::log(domainVector[i] - min_val) -
           std::log(max_val - domainVector[i]);
     }
     else if (queso_isfinite(min_val) &&
              !queso_isfinite(max_val)) {
 
       if (domainVector[i] == min_val) {
         // Exit early if we can
         return -INFINITY;
       }
 
       // Left-hand side finite, but right-hand side is not.
       // Do only left-hand transform.
       transformedDomainVector[i] = std::log(domainVector[i] - min_val);
 
       lnjacobian += -std::log(domainVector[i] - min_val);
     }
     else if (!queso_isfinite(min_val) &&
              queso_isfinite(max_val)) {
 
       if (domainVector[i] == max_val) {
         // Exit early if we can
         return -INFINITY;
       }
 
       // Right-hand side is finite, but left-hand side is not.
       // Do only right-hand transform.
       transformedDomainVector[i] = -std::log(max_val - domainVector[i]);
 
       lnjacobian += -std::log(max_val - domainVector[i]);
     }
     else {
       // No transform.
       transformedDomainVector[i] = domainVector[i];
     }
   }
 
   V diffVec(transformedDomainVector - this->lawExpVector());
   if (m_diagonalCovMatrix) {
     returnValue = ((diffVec * diffVec) /
         this->lawVarVector()).sumOfComponents();
     if (m_normalizationStyle == 0) {
       unsigned int iMax = this->lawVarVector().sizeLocal();
       for (unsigned int i = 0; i < iMax; ++i) {
         lnDeterminant += log(this->lawVarVector()[i]);
       }
     }
   }
   else {
     V tmpVec = this->m_lawCovMatrix->invertMultiply(diffVec);
     returnValue = (diffVec * tmpVec).sumOfComponents();
     if (m_normalizationStyle == 0) {
       lnDeterminant = this->m_lawCovMatrix->lnDeterminant();
     }
   }
   if (m_normalizationStyle == 0) {
     returnValue += ((double) this->lawVarVector().sizeLocal()) * log(2 * M_PI);
     returnValue += lnDeterminant;
   }
   returnValue *= -0.5;
   returnValue += m_logOfNormalizationFactor;
   returnValue += lnjacobian;
 
   return returnValue;
 }

template<class V , class M >

void QUESO::InvLogitGaussianJointPdf< V, M >::print ( std::ostream & os ) const

virtual

Prints the distribution.

When the mean and variance are printed, what's printed is the mean and variance of the underlying Gaussian distribution, not the output of the distributionMean and distributionVariance methods.

Reimplemented from QUESO::BaseJointPdf< V, M >.

Definition at line 91 of file InvLogitGaussianJointPdf.C.

 {
   // Print m_env?
 
   os << "Start printing InvLogitGaussianJointPdf<V, M>" << std::endl;
   os << "m_prefix:" << std::endl;
   os << this->m_prefix << std::endl;
   os << "m_domainSet:" << std::endl;
   os << this->m_domainBoxSubset << std::endl;
   os << "m_normalizationStyle:" << std::endl;
   os << this->m_normalizationStyle << std::endl;
   os << "m_logOfNormalizationFactor:" << std::endl;
   os << this->m_logOfNormalizationFactor << std::endl;
   os << "Mean:" << std::endl;
   os << this->lawExpVector() << std::endl;
   os << "Variance vector:" << std::endl;
   os << this->lawVarVector() << std::endl;
   os << "Covariance matrix:" << std::endl;
   os << this->lawCovMatrix() << std::endl;
   os << "Diagonal covariance?" << std::endl;
   os << this->m_diagonalCovMatrix << std::endl;
   os << "End printing InvLogitGaussianJointPdf<V, M>" << std::endl;
 }

template<class V , class M >

void QUESO::InvLogitGaussianJointPdf< V, M >::updateLawCovMatrix ( const M & newLawCovMatrix )

Updates the lower triangular matrix from Cholesky decomposition of the covariance matrix to the new value newLowerCholLawCovMatrix.

This method deletes old expected values (allocated at construction or last call to this method).

Definition at line 286 of file InvLogitGaussianJointPdf.C.

 {
   // delete old expected values (allocated at construction or last call to this function)
   delete m_lawCovMatrix;
   m_lawCovMatrix = new M(newLawCovMatrix);
 }

template<class V , class M >

void QUESO::InvLogitGaussianJointPdf< V, M >::updateLawExpVector ( const V & newLawExpVector )

Updates the mean of the Gaussian (not transformed) with the new value newLawExpVector.

This method deletes old expected values (allocated at construction or last call to this method).

Definition at line 277 of file InvLogitGaussianJointPdf.C.

 {
   // delete old expected values (allocated at construction or last call to this function)
   delete m_lawExpVector;
   m_lawExpVector = new V(newLawExpVector);
 }