25 #include <queso/LogNormalJointPdf.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)+
"gau").c_str(),domainSet),
40 m_lawExpVector (new V(lawExpVector)),
41 m_lawVarVector (new V(lawVarVector)),
42 m_diagonalCovMatrix(true)
65 template<
class V,
class M>
68 delete m_lawVarVector;
69 delete m_lawExpVector;
72 template <
class V,
class M>
76 return *m_lawExpVector;
79 template <
class V,
class M>
83 return *m_lawVarVector;
86 template<
class V,
class M>
89 const V& domainVector,
90 const V* domainDirection,
93 V* hessianEffect)
const
95 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
96 *m_env.subDisplayFile() <<
"Entering LogNormalJointPdf<V,M>::actualValue()"
97 <<
", meanVector = " << *m_lawExpVector
98 <<
": domainVector = " << domainVector
99 <<
", domainVector.sizeLocal() = " << domainVector.sizeLocal()
100 <<
", this->m_domainSet.vectorSpace().dimLocal() = " << this->m_domainSet.vectorSpace().dimLocal()
106 queso_require_msg(!(hessianMatrix || hessianEffect),
"incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
108 double returnValue = 0.;
110 V zeroVector(domainVector);
111 zeroVector.cwSet(0.);
112 if (domainVector.atLeastOneComponentSmallerOrEqualThan(zeroVector)) {
116 else if (this->m_domainSet.contains(domainVector) ==
false) {
122 returnValue = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
125 (*gradVector) *= returnValue;
129 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
130 *m_env.subDisplayFile() <<
"Leaving LogNormalJointPdf<V,M>::actualValue()"
131 <<
", meanVector = " << *m_lawExpVector
132 <<
": domainVector = " << domainVector
133 <<
", returnValue = " << returnValue
140 template<
class V,
class M>
143 const V& domainVector,
144 const V* domainDirection,
147 V* hessianEffect)
const
149 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
150 *m_env.subDisplayFile() <<
"Entering LogNormalJointPdf<V,M>::lnValue()"
151 <<
", meanVector = " << *m_lawExpVector
152 <<
": domainVector = " << domainVector
156 queso_require_msg(!(domainDirection || hessianMatrix || hessianEffect),
"incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
158 double returnValue = 0.;
160 V zeroVector(domainVector);
161 zeroVector.cwSet(0.);
162 if (domainVector.atLeastOneComponentSmallerOrEqualThan(zeroVector)) {
164 returnValue = -INFINITY;
166 else if (this->m_domainSet.contains(domainVector) ==
false) {
168 returnValue = -INFINITY;
171 if (m_diagonalCovMatrix) {
172 V diffVec(zeroVector);
173 for (
unsigned int i = 0; i < domainVector.sizeLocal(); ++i) {
174 diffVec[i] = std::log(domainVector[i]) - this->lawExpVector()[i];
182 (*gradVector)[i] = -(1.0 / domainVector[i]) -
183 diffVec[i] / (domainVector[i] * this->lawVarVector()[i]);
186 returnValue = ((diffVec*diffVec)/this->lawVarVector()).sumOfComponents();
189 for (
unsigned int i = 0; i < domainVector.sizeLocal(); ++i) {
190 returnValue -= std::log(domainVector[i]);
191 if (m_normalizationStyle == 0) {
192 returnValue -= std::log(std::sqrt(2. * M_PI * this->lawVarVector()[i]));
197 queso_error_msg(
"situation with a non-diagonal covariance matrix makes no sense");
199 returnValue += m_logOfNormalizationFactor;
202 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
203 *m_env.subDisplayFile() <<
"Leaving LogNormalJointPdf<V,M>::lnValue()"
204 <<
", meanVector = " << *m_lawExpVector
205 <<
": domainVector = " << domainVector
206 <<
", returnValue = " << returnValue
213 template<
class V,
class M>
220 if (m_diagonalCovMatrix) {
221 unsigned int n_params = meanVector.sizeLocal();
224 for (
unsigned int i = 0; i < n_params; ++i) {
225 meanVector[i] = std::exp(this->lawExpVector()[i] + 0.5*this->lawVarVector()[i]);
229 queso_error_msg(
"situation with a non-diagonal covariance matrix makes no sense");
233 template<
class V,
class M>
240 if (m_diagonalCovMatrix) {
241 unsigned int n_params = this->lawExpVector().sizeLocal();
246 covMatrix.zeroLower();
247 covMatrix.zeroUpper();
249 for (
unsigned int i = 0; i < n_params; ++i) {
250 covMatrix(i,i) = (std::exp(this->lawVarVector()[i]) - 1) *
251 std::exp(2*this->lawExpVector()[i] + this->lawVarVector()[i]);
255 queso_error_msg(
"situation with a non-diagonal covariance matrix makes no sense");
259 template<
class V,
class M>
265 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
266 *m_env.subDisplayFile() <<
"Entering LogNormalJointPdf<V,M>::computeLogOfNormalizationFactor()"
270 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
271 *m_env.subDisplayFile() <<
"Leaving LogNormalJointPdf<V,M>::computeLogOfNormalizationFactor()"
272 <<
", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
~LogNormalJointPdf()
Destructor.
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the Log-Normal PDF (scalar function).
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
double lnValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Logarithm of the value of the Log-Normal PDF (scalar function).
virtual void distributionVariance(M &covMatrix) const
Covariance matrix of the underlying random variable.
#define queso_require_equal_to_msg(expr1, expr2, msg)
const V & lawExpVector() const
Access to the vector of mean values and private attribute: m_lawExpVector.
A templated (base) class for handling joint PDFs.
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
#define queso_assert_equal_to(expr1, expr2)
A templated class for handling sets.
A class for handling Log-Normal joint PDFs.
double commonComputeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Common method (to the derived classes) to compute the logarithm of the normalization factor...
unsigned int displayVerbosity() const
virtual void distributionMean(V &meanVector) const
Mean value of the underlying random variable.
#define queso_require_msg(asserted, msg)
const BaseEnvironment & m_env
LogNormalJointPdf(const char *prefix, const VectorSet< V, M > &domainSet, const V &lawExpVector, const V &lawVarVector)
Constructor.
const V & lawVarVector() const
Access to the vector of variance values and private attribute: m_lawVarVector.
#define queso_error_msg(msg)