25 #include <queso/InvLogitGaussianJointPdf.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
32 template<
class V,
class M>
36 const V& lawExpVector,
37 const V& lawVarVector)
39 BaseJointPdf<V,M>(((std::string)(prefix)+
"invlogit_gau").c_str(),
41 m_lawExpVector(new V(lawExpVector)),
42 m_lawVarVector(new V(lawVarVector)),
43 m_diagonalCovMatrix(true),
44 m_lawCovMatrix(m_domainSet.vectorSpace().newDiagMatrix(lawVarVector)),
45 m_domainBoxSubset(domainBoxSubset)
50 template<
class V,
class M>
54 const V& lawExpVector,
55 const M& lawCovMatrix)
57 BaseJointPdf<V,M>(((std::string)(prefix)+
"invlogit_gau").c_str(),
59 m_lawExpVector(new V(lawExpVector)),
60 m_lawVarVector(domainBoxSubset.vectorSpace().newVector(INFINITY)),
61 m_diagonalCovMatrix(false),
62 m_lawCovMatrix(new M(lawCovMatrix)),
63 m_domainBoxSubset(domainBoxSubset)
67 template<
class V,
class M>
70 delete m_lawCovMatrix;
71 delete m_lawVarVector;
72 delete m_lawExpVector;
75 template <
class V,
class M>
79 return *m_lawExpVector;
82 template <
class V,
class M>
86 return *m_lawVarVector;
89 template<
class V,
class M>
92 const V& domainVector,
93 const V* domainDirection,
96 V* hessianEffect)
const
100 returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
105 template<
class V,
class M>
108 const V& domainVector,
109 const V* domainDirection,
112 V* hessianEffect)
const
115 double lnDeterminant = 0.0;
116 V transformedDomainVector(domainVector);
118 V min_domain_bounds(this->m_domainBoxSubset.minValues());
119 V max_domain_bounds(this->m_domainBoxSubset.maxValues());
121 double lnjacobian = 0.0;
122 for (
unsigned int i = 0; i < domainVector.sizeLocal(); i++) {
123 double min_val = min_domain_bounds[i];
124 double max_val = max_domain_bounds[i];
126 if (boost::math::isfinite(min_val) &&
127 boost::math::isfinite(max_val)) {
129 if (domainVector[i] == min_val || domainVector[i] == max_val) {
135 transformedDomainVector[i] = std::log(domainVector[i] - min_val) -
136 std::log(max_val - domainVector[i]);
138 lnjacobian += std::log(max_val - min_val) -
139 std::log(domainVector[i] - min_val) -
140 std::log(max_val - domainVector[i]);
142 else if (boost::math::isfinite(min_val) &&
143 !boost::math::isfinite(max_val)) {
145 if (domainVector[i] == min_val) {
152 transformedDomainVector[i] = std::log(domainVector[i] - min_val);
154 lnjacobian += -std::log(domainVector[i] - min_val);
156 else if (!boost::math::isfinite(min_val) &&
157 boost::math::isfinite(max_val)) {
159 if (domainVector[i] == max_val) {
166 transformedDomainVector[i] = -std::log(max_val - domainVector[i]);
168 lnjacobian += -std::log(max_val - domainVector[i]);
172 transformedDomainVector[i] = domainVector[i];
176 V diffVec(transformedDomainVector - this->lawExpVector());
177 if (m_diagonalCovMatrix) {
178 returnValue = ((diffVec * diffVec) /
179 this->lawVarVector()).sumOfComponents();
180 if (m_normalizationStyle == 0) {
181 unsigned int iMax = this->lawVarVector().sizeLocal();
182 for (
unsigned int i = 0; i < iMax; ++i) {
183 lnDeterminant += log(this->lawVarVector()[i]);
188 V tmpVec = this->m_lawCovMatrix->invertMultiply(diffVec);
189 returnValue = (diffVec * tmpVec).sumOfComponents();
190 if (m_normalizationStyle == 0) {
191 lnDeterminant = this->m_lawCovMatrix->lnDeterminant();
194 if (m_normalizationStyle == 0) {
195 returnValue += ((double) this->lawVarVector().sizeLocal()) * log(2 * M_PI);
196 returnValue += lnDeterminant;
199 returnValue += m_logOfNormalizationFactor;
200 returnValue += lnjacobian;
205 template<
class V,
class M>
211 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
212 *m_env.subDisplayFile() <<
"Entering GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
216 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
217 *m_env.subDisplayFile() <<
"Leaving GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
218 <<
", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
225 template<
class V,
class M>
230 delete m_lawExpVector;
231 m_lawExpVector =
new V(newLawExpVector);
234 template<
class V,
class M>
239 delete m_lawCovMatrix;
240 m_lawCovMatrix =
new M(newLawCovMatrix);
243 template<
class V,
class M>
247 return *m_lawCovMatrix;
void updateLawCovMatrix(const M &newLawCovMatrix)
Updates the lower triangular matrix from Cholesky decomposition of the covariance matrix to the new v...
const V & lawVarVector() const
Access to the vector of variance values and private attribute: m_lawVarVector.
~InvLogitGaussianJointPdf()
Destructor.
InvLogitGaussianJointPdf(const char *prefix, const BoxSubset< V, M > &domainBoxSubset, const V &lawExpVector, const V &lawVarVector)
Constructor.
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).
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the (transformed) Gaussian PDF.
Class representing a subset of a vector space shaped like a hypercube.
double commonComputeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Common method (to the derived classes) to compute the logarithm of the normalization factor...
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
const M & lawCovMatrix() const
Returns the covariance matrix; access to protected attribute m_lawCovMatrix.
A class for handling hybrid (transformed) Gaussians with bounds.
const V & lawExpVector() const
Access to the vector of mean values of the Gaussian (not transformed) and private attribute: m_lawExp...
A templated (base) class for handling joint PDFs.
void updateLawExpVector(const V &newLawExpVector)
Updates the mean of the Gaussian (not transformed) with the new value newLawExpVector.