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 if (m_normalizationStyle == 0) {
190 for (
unsigned int i = 0; i < domainVector.sizeLocal(); ++i) {
191 returnValue -= std::log(domainVector[i] * std::sqrt(2. * M_PI * this->lawVarVector()[i]));
196 queso_error_msg(
"situation with a non-diagonal covariance matrix makes no sense");
198 returnValue += m_logOfNormalizationFactor;
201 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 55)) {
202 *m_env.subDisplayFile() <<
"Leaving LogNormalJointPdf<V,M>::lnValue()"
203 <<
", meanVector = " << *m_lawExpVector
204 <<
": domainVector = " << domainVector
205 <<
", returnValue = " << returnValue
212 template<
class V,
class M>
219 if (m_diagonalCovMatrix) {
220 unsigned int n_params = meanVector.sizeLocal();
223 for (
unsigned int i = 0; i < n_params; ++i) {
224 meanVector[i] = std::exp(this->lawExpVector()[i] + 0.5*this->lawVarVector()[i]);
228 queso_error_msg(
"situation with a non-diagonal covariance matrix makes no sense");
232 template<
class V,
class M>
239 if (m_diagonalCovMatrix) {
240 unsigned int n_params = this->lawExpVector().sizeLocal();
245 covMatrix.zeroLower();
246 covMatrix.zeroUpper();
248 for (
unsigned int i = 0; i < n_params; ++i) {
249 covMatrix(i,i) = (std::exp(this->lawVarVector()[i]) - 1) *
250 std::exp(2*this->lawExpVector()[i] + this->lawVarVector()[i]);
254 queso_error_msg(
"situation with a non-diagonal covariance matrix makes no sense");
258 template<
class V,
class M>
264 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
265 *m_env.subDisplayFile() <<
"Entering LogNormalJointPdf<V,M>::computeLogOfNormalizationFactor()"
269 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
270 *m_env.subDisplayFile() <<
"Leaving LogNormalJointPdf<V,M>::computeLogOfNormalizationFactor()"
271 <<
", 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)