A class for handling Gaussian joint PDFs. More...

#include <GaussianJointPdf.h>

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

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

Public Member Functions
virtual void	print (std::ostream &os) const
	Print method for informational and logging purposes. More...

Constructor/Destructor methods
	GaussianJointPdf (const char *prefix, const VectorSet< V, M > &domainSet, const V &lawExpVector, const V &lawVarVector)
	Constructor. More...

	GaussianJointPdf (const char *prefix, const VectorSet< V, M > &domainSet, const V &lawExpVector, const M &lawCovMatrix)
	Constructor. More...

	~GaussianJointPdf ()
	Destructor. More...

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

double	lnValue (const V &domainVector, const V domainDirection, V gradVector, M hessianMatrix, V hessianEffect) const
	Logarithm of the value of the 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 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 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...

Protected Attributes
V *	m_lawExpVector

V *	m_lawVarVector

bool	m_diagonalCovMatrix

const M *	m_lawCovMatrix

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

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

Detailed Description

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

A class for handling Gaussian joint PDFs.

This class allows the mathematical definition of a Gaussian Joint PDF.

Definition at line 50 of file GaussianJointPdf.h.

Constructor & Destructor Documentation

template<class V , class M >

QUESO::GaussianJointPdf< V, M >::GaussianJointPdf	(	const char *	prefix,
		const VectorSet< V, M > &	domainSet,
		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, and a vector of covariance values lawVarVector (an alternative representation for a diagonal covariance matrix).

Definition at line 33 of file GaussianJointPdf.C.

References QUESO::BaseEnvironment::displayVerbosity(), QUESO::GaussianJointPdf< V, M >::lawExpVector(), QUESO::GaussianJointPdf< V, M >::lawVarVector(), QUESO::BaseScalarFunction< V, M >::m_env, QUESO::BaseScalarFunction< V, M >::m_prefix, and QUESO::BaseEnvironment::subDisplayFile().

   :
   BaseJointPdf<V,M>(((std::string)(prefix)+"gau").c_str(),domainSet),
   m_lawExpVector     (new V(lawExpVector)),
   m_lawVarVector     (new V(lawVarVector)),
   m_diagonalCovMatrix(true),
   m_lawCovMatrix     (m_domainSet.vectorSpace().newDiagMatrix(lawVarVector))
 {
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
     *m_env.subDisplayFile() << "Entering GaussianJointPdf<V,M>::constructor() [1]"
                             << ": prefix = " << m_prefix
                             << std::endl;
   }
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
     *m_env.subDisplayFile() << "In GaussianJointPdf<V,M>::constructor()"
                           //<< ", prefix = "     << m_prefix
                             << ": meanVector = " << this->lawExpVector()
                       << ", Variances = "  << this->lawVarVector()
                             << std::endl;
   }
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
     *m_env.subDisplayFile() << "Leaving GaussianJointPdf<V,M>::constructor() [1]"
                             << ": prefix = " << m_prefix
                             << std::endl;
   }
 }

template<class V , class M >

QUESO::GaussianJointPdf< V, M >::GaussianJointPdf	(	const char *	prefix,
		const VectorSet< V, M > &	domainSet,
		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, and a covariance matrix, lawCovMatrix.

Definition at line 68 of file GaussianJointPdf.C.

References QUESO::BaseEnvironment::displayVerbosity(), QUESO::GaussianJointPdf< V, M >::lawExpVector(), QUESO::BaseScalarFunction< V, M >::m_env, QUESO::BaseScalarFunction< V, M >::m_prefix, and QUESO::BaseEnvironment::subDisplayFile().

   :
   BaseJointPdf<V,M>(((std::string)(prefix)+"gau").c_str(),domainSet),
   m_lawExpVector     (new V(lawExpVector)),
   m_lawVarVector     (domainSet.vectorSpace().newVector(INFINITY)), // FIX ME
   m_diagonalCovMatrix(false),
   m_lawCovMatrix     (new M(lawCovMatrix))
 {
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
     *m_env.subDisplayFile() << "Entering GaussianJointPdf<V,M>::constructor() [2]"
                             << ": prefix = " << m_prefix
                             << std::endl;
   }
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
     *m_env.subDisplayFile() << "In GaussianJointPdf<V,M>::constructor()"
                           //<< ", prefix = "            << m_prefix
                             << ": meanVector = "        << this->lawExpVector()
                       << ", Covariance Matrix = " << lawCovMatrix
                             << std::endl;
   }
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
     *m_env.subDisplayFile() << "Leaving GaussianJointPdf<V,M>::constructor() [2]"
                             << ": prefix = " << m_prefix
                             << std::endl;
   }
 }

template<class V , class M >

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

Destructor.

Definition at line 102 of file GaussianJointPdf.C.

 {
   delete m_lawCovMatrix;
   delete m_lawVarVector;
   delete m_lawExpVector;
 }

Member Function Documentation

template<class V , class M >

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

virtual

Actual value of the Gaussian PDF:

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

Implements QUESO::BaseJointPdf< V, M >.

Definition at line 152 of file GaussianJointPdf.C.

References queso_require_equal_to_msg, and queso_require_msg.

 {
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
     *m_env.subDisplayFile() << "Entering GaussianJointPdf<V,M>::actualValue()"
                             << ", meanVector = "   << *m_lawExpVector
                       << ", lawCovMatrix = " << *m_lawCovMatrix
                             << ": domainVector = " << domainVector
                             << std::endl;
   }
 
   queso_require_equal_to_msg(domainVector.sizeLocal(), this->m_domainSet.vectorSpace().dimLocal(), "invalid input");
 
   queso_require_msg(!(hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
 
   double returnValue = 0.;
 
   if (this->m_domainSet.contains(domainVector) == false) {
     // What should the gradient be here?
     returnValue = 0.;
   }
   else {
     // Already normalised (so the gradient will have the normalisation constant
     // in it)
     returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
 
     if (gradVector) {
       (*gradVector) *= returnValue;
     }
   }
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
     *m_env.subDisplayFile() << "Leaving GaussianJointPdf<V,M>::actualValue()"
                             << ", meanVector = "   << *m_lawExpVector
                       << ", lawCovMatrix = " << *m_lawCovMatrix
                             << ": domainVector = " << domainVector
                             << ", returnValue = "  << returnValue
                             << std::endl;
   }
 
   return returnValue;
 }

template<class V , class M >

double QUESO::GaussianJointPdf< 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 329 of file GaussianJointPdf.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::GaussianJointPdf< V, M >::distributionMean ( V & meanVector ) const

virtual

Mean value of the underlying random variable.

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

Definition at line 301 of file GaussianJointPdf.C.

 {
   meanVector = this->lawExpVector();
 }

template<class V , class M >

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

virtual

Covariance matrix of the underlying random variable.

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

Definition at line 308 of file GaussianJointPdf.C.

References queso_assert_equal_to.

 {
   queso_assert_equal_to (covMatrix.numCols(), covMatrix.numRowsGlobal());
 
   if (m_diagonalCovMatrix) {
     covMatrix.zeroLower();
     covMatrix.zeroUpper();
 
     unsigned int n_comp = this->lawVarVector().sizeLocal();
     queso_assert_equal_to (n_comp, covMatrix.numCols());
 
     for (unsigned int i = 0; i < n_comp; ++i) {
       covMatrix(i,i) = this->lawVarVector()[i];
     }
   } else {
     covMatrix = *this->m_lawCovMatrix;
   }
 }

template<class V , class M >

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

Returns the covariance matrix; access to protected attribute m_lawCovMatrix.

Definition at line 369 of file GaussianJointPdf.C.

Referenced by QUESO::ScaledCovMatrixTKGroup< V, M >::setPreComputingPosition(), and QUESO::MetropolisAdjustedLangevinTK< V, M >::setPreComputingPosition().

 {
   return *m_lawCovMatrix;
 }

template<class V , class M >

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

Access to the vector of mean values and private attribute: m_lawExpVector.

Definition at line 111 of file GaussianJointPdf.C.

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

 {
   return *m_lawExpVector;
 }

template<class V , class M >

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

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

Definition at line 118 of file GaussianJointPdf.C.

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

 {
   return *m_lawVarVector;
 }

template<class V , class M >

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

virtual

Logarithm of the value of the 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},$

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 201 of file GaussianJointPdf.C.

References queso_require_msg.

 {
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
     *m_env.subDisplayFile() << "Entering GaussianJointPdf<V,M>::lnValue()"
                             << ", meanVector = "   << *m_lawExpVector
                       << ", lawCovMatrix = " << *m_lawCovMatrix
                             << ": domainVector = " << domainVector
                             << std::endl;
   }
 
   queso_require_msg(!(domainDirection || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
 
   double returnValue = 0.;
 
   double lnDeterminant = 0.;
   if (this->m_domainSet.contains(domainVector) == false) {
     // What should the gradient be here?
     returnValue = -INFINITY;
   }
   else {
     V diffVec(domainVector - this->lawExpVector());
     if (m_diagonalCovMatrix) {
       returnValue = ((diffVec*diffVec)/this->lawVarVector()).sumOfComponents();
 
       // Compute the gradient of log of the pdf.
       // The log of a Gaussian pdf is:
       // f(x) = - 1/2 (x - \mu)^T \Sigma^{-1} (x - \mu)
       // Therefore
       // \frac{df}{dx}(x) = - (x - \mu)^T \Sigma^{-1}
       //                  = - (\Sigma^{-1}^T (x - \mu))^T
       //                  = - (\Sigma^{-1} (x - \mu))^T
       //                  = - \Sigma^{-1} (x - \mu)  (row/column vector doesn't matter)
       //
       // So if \Sigma is diagonal we have a component-wise product of two
       // vectors (x - \mu) and the diagonal elements of \Sigma^{-1}
       if (gradVector) {
         (*gradVector) = diffVec;  // Copy
         (*gradVector) /= this->lawVarVector();
         (*gradVector) *= -1.0;
       }
 
       if (m_normalizationStyle == 0) {
         unsigned int iMax = this->lawVarVector().sizeLocal();
         for (unsigned int i = 0; i < iMax; ++i) {
           lnDeterminant += std::log(this->lawVarVector()[i]);
         }
       }
     }
     else {
       V tmpVec = this->m_lawCovMatrix->invertMultiply(diffVec);
       returnValue = (diffVec*tmpVec).sumOfComponents();
 
       // Compute the gradient of log of the pdf.
       // The log of a Gaussian pdf is:
       // f(x) = - 1/2 (x - \mu)^T \Sigma^{-1} (x - \mu)
       // Therefore
       // \frac{df}{dx}(x) = - (x - \mu)^T \Sigma^{-1}
       //                  = - (\Sigma^{-1}^T (x - \mu))^T
       //                  = - (\Sigma^{-1} (x - \mu))^T
       //                  = - \Sigma^{-1} (x - \mu)  (row/column vector doesn't matter)
       if (gradVector) {
         (*gradVector) = tmpVec;
         (*gradVector) *= -1.0;
       }
 
       if (m_normalizationStyle == 0) {
         lnDeterminant = this->m_lawCovMatrix->lnDeterminant();
       }
     }
     if (m_normalizationStyle == 0) {
       returnValue += ((double) this->lawVarVector().sizeLocal()) * std::log(2*M_PI);   // normalization of pdf
       returnValue += lnDeterminant; // normalization of pdf
     }
     returnValue *= -0.5;
   }
   returnValue += m_logOfNormalizationFactor;
 
   if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
     *m_env.subDisplayFile() << "Leaving GaussianJointPdf<V,M>::lnValue()"
                             << ", m_normalizationStyle = " << m_normalizationStyle
                             << ", m_diagonalCovMatrix = " << m_diagonalCovMatrix
                             << ", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
                             << ", lnDeterminant = " << lnDeterminant
                             << ", meanVector = "           << *m_lawExpVector
                             << ", lawCovMatrix = "         << *m_lawCovMatrix
                             << ": domainVector = "         << domainVector
                             << ", returnValue = "          << returnValue
                             << std::endl;
   }
 
   return returnValue;
 }

template<class V , class M >

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

virtual

Print method for informational and logging purposes.

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

Definition at line 125 of file GaussianJointPdf.C.

 {
   // Print m_env?
 
   os << "Start printing GaussianJointPdf<V, M>" << std::endl;
   os << "m_prefix:" << std::endl;
   os << this->m_prefix << std::endl;
   os << "m_domainSet:" << std::endl;
   os << this->m_domainSet << 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 GaussianJointPdf<V, M>" << std::endl;
 }

template<class V , class M >

void QUESO::GaussianJointPdf< 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 359 of file GaussianJointPdf.C.

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

template<class V , class M >

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

Updates the mean with the new value newLawExpVector.

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

Definition at line 349 of file GaussianJointPdf.C.

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

Member Data Documentation

template<class V = GslVector, class M = GslMatrix>

bool QUESO::GaussianJointPdf< V, M >::m_diagonalCovMatrix

protected

Definition at line 126 of file GaussianJointPdf.h.

template<class V = GslVector, class M = GslMatrix>

const M* QUESO::GaussianJointPdf< V, M >::m_lawCovMatrix

protected

Definition at line 127 of file GaussianJointPdf.h.

template<class V = GslVector, class M = GslMatrix>

V* QUESO::GaussianJointPdf< V, M >::m_lawExpVector

protected

Definition at line 124 of file GaussianJointPdf.h.

template<class V = GslVector, class M = GslMatrix>

V* QUESO::GaussianJointPdf< V, M >::m_lawVarVector

protected

Definition at line 125 of file GaussianJointPdf.h.

The documentation for this class was generated from the following files:

src/stats/inc/GaussianJointPdf.h
src/stats/src/GaussianJointPdf.C

Public Member Functions

Protected Attributes

Additional Inherited Members

Detailed Description

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

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation

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