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
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...

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 91 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 207 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 >

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

Returns the covariance matrix; access to protected attribute m_lawCovMatrix.

Definition at line 245 of file InvLogitGaussianJointPdf.C.

Referenced by QUESO::MetropolisHastingsSG< P_V, P_M >::alpha(), and 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.

Referenced by QUESO::MetropolisHastingsSG< P_V, P_M >::alpha().

 {
   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.

Referenced by QUESO::MetropolisHastingsSG< P_V, P_M >::alpha().

 {
   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 107 of file InvLogitGaussianJointPdf.C.

 {
   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 (boost::math::isfinite(min_val) &&
         boost::math::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 (boost::math::isfinite(min_val) &&
              !boost::math::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 (!boost::math::isfinite(min_val) &&
              boost::math::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 >::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 236 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 227 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);
 }