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>
95 os <<
"Start printing InvLogitGaussianJointPdf<V, M>" << std::endl;
96 os <<
"m_prefix:" << std::endl;
97 os << this->m_prefix << std::endl;
98 os <<
"m_domainSet:" << std::endl;
99 os << this->m_domainBoxSubset << std::endl;
100 os <<
"m_normalizationStyle:" << std::endl;
101 os << this->m_normalizationStyle << std::endl;
102 os <<
"m_logOfNormalizationFactor:" << std::endl;
103 os << this->m_logOfNormalizationFactor << std::endl;
104 os <<
"Mean:" << std::endl;
105 os << this->lawExpVector() << std::endl;
106 os <<
"Variance vector:" << std::endl;
107 os << this->lawVarVector() << std::endl;
108 os <<
"Covariance matrix:" << std::endl;
109 os << this->lawCovMatrix() << std::endl;
110 os <<
"Diagonal covariance?" << std::endl;
111 os << this->m_diagonalCovMatrix << std::endl;
112 os <<
"End printing InvLogitGaussianJointPdf<V, M>" << std::endl;
115 template<
class V,
class M>
118 const V& domainVector,
119 const V* domainDirection,
122 V* hessianEffect)
const
126 returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
131 template<
class V,
class M>
134 const V& domainVector,
135 const V* domainDirection,
138 V* hessianEffect)
const
141 double lnDeterminant = 0.0;
142 V transformedDomainVector(domainVector);
144 V min_domain_bounds(this->m_domainBoxSubset.minValues());
145 V max_domain_bounds(this->m_domainBoxSubset.maxValues());
147 double lnjacobian = 0.0;
148 for (
unsigned int i = 0; i < domainVector.sizeLocal(); i++) {
149 double min_val = min_domain_bounds[i];
150 double max_val = max_domain_bounds[i];
155 if (domainVector[i] == min_val || domainVector[i] == max_val) {
161 transformedDomainVector[i] = std::log(domainVector[i] - min_val) -
162 std::log(max_val - domainVector[i]);
164 lnjacobian += std::log(max_val - min_val) -
165 std::log(domainVector[i] - min_val) -
166 std::log(max_val - domainVector[i]);
171 if (domainVector[i] == min_val) {
178 transformedDomainVector[i] = std::log(domainVector[i] - min_val);
180 lnjacobian += -std::log(domainVector[i] - min_val);
185 if (domainVector[i] == max_val) {
192 transformedDomainVector[i] = -std::log(max_val - domainVector[i]);
194 lnjacobian += -std::log(max_val - domainVector[i]);
198 transformedDomainVector[i] = domainVector[i];
202 V diffVec(transformedDomainVector - this->lawExpVector());
203 if (m_diagonalCovMatrix) {
204 returnValue = ((diffVec * diffVec) /
205 this->lawVarVector()).sumOfComponents();
206 if (m_normalizationStyle == 0) {
207 unsigned int iMax = this->lawVarVector().sizeLocal();
208 for (
unsigned int i = 0; i < iMax; ++i) {
209 lnDeterminant += log(this->lawVarVector()[i]);
214 V tmpVec = this->m_lawCovMatrix->invertMultiply(diffVec);
215 returnValue = (diffVec * tmpVec).sumOfComponents();
216 if (m_normalizationStyle == 0) {
217 lnDeterminant = this->m_lawCovMatrix->lnDeterminant();
220 if (m_normalizationStyle == 0) {
221 returnValue += ((double) this->lawVarVector().sizeLocal()) * log(2 * M_PI);
222 returnValue += lnDeterminant;
225 returnValue += m_logOfNormalizationFactor;
226 returnValue += lnjacobian;
231 template<
class V,
class M>
239 queso_not_implemented();
243 template<
class V,
class M>
251 queso_not_implemented();
255 template<
class V,
class M>
261 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
262 *m_env.subDisplayFile() <<
"Entering GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
266 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
267 *m_env.subDisplayFile() <<
"Leaving GaussianJointPdf<V,M>::computeLogOfNormalizationFactor()"
268 <<
", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
275 template<
class V,
class M>
280 delete m_lawExpVector;
281 m_lawExpVector =
new V(newLawExpVector);
284 template<
class V,
class M>
289 delete m_lawCovMatrix;
290 m_lawCovMatrix =
new M(newLawCovMatrix);
293 template<
class V,
class M>
297 return *m_lawCovMatrix;
InvLogitGaussianJointPdf(const char *prefix, const BoxSubset< V, M > &domainBoxSubset, const V &lawExpVector, const V &lawVarVector)
Constructor.
Class representing a subset of a vector space shaped like a hypercube.
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the (transformed) Gaussian PDF.
void updateLawExpVector(const V &newLawExpVector)
Updates the mean of the Gaussian (not transformed) with the new value newLawExpVector.
const M & lawCovMatrix() const
Returns the covariance matrix; access to protected attribute m_lawCovMatrix.
A class for handling hybrid (transformed) Gaussians with bounds.
A templated (base) class for handling joint PDFs.
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
~InvLogitGaussianJointPdf()
Destructor.
const V & lawExpVector() const
Access to the vector of mean values of the Gaussian (not transformed) and private attribute: m_lawExp...
double commonComputeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Common method (to the derived classes) to compute the logarithm of the normalization factor...
virtual void print(std::ostream &os) const
Prints the distribution.
void updateLawCovMatrix(const M &newLawCovMatrix)
Updates the lower triangular matrix from Cholesky decomposition of the covariance matrix to the new v...
virtual void distributionVariance(M &covMatrix) const
Covariance matrix of the underlying random variable.
virtual void distributionMean(V &meanVector) const
Mean value of the underlying random variable.
const V & lawVarVector() const
Access to the vector of variance values and private attribute: m_lawVarVector.
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).
bool queso_isfinite(T arg)