25 #include <queso/BayesianJointPdf.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
32 template<
class V,
class M>
37 double likelihoodExponent,
40 BaseJointPdf<V,M>(((std::string)(prefix)+
"bay").c_str(),intersectionDomain),
41 m_priorDensity (priorDensity),
42 m_likelihoodFunction (likelihoodFunction),
43 m_likelihoodExponent (likelihoodExponent),
44 m_lastComputedLogPrior (0.),
45 m_lastComputedLogLikelihood(0.),
46 m_tmpVector1 (m_domainSet.vectorSpace().zeroVector()),
47 m_tmpVector2 (m_domainSet.vectorSpace().zeroVector()),
48 m_tmpMatrix (m_domainSet.vectorSpace().newMatrix())
52 template<
class V,
class M>
58 template<
class V,
class M>
62 m_priorDensity.setNormalizationStyle(value);
66 template<
class V,
class M>
70 return m_lastComputedLogPrior;
73 template<
class V,
class M>
77 return m_lastComputedLogLikelihood;
80 template<
class V,
class M>
83 const V& domainVector,
84 const V* domainDirection,
87 V* hessianEffect)
const
89 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
90 *m_env.subDisplayFile() <<
"Entering BayesianJointPdf<V,M>::actualValue()"
91 <<
": domainVector = " << domainVector
97 "BayesianJointPdf<V,M>::actualValue()",
101 if (gradVector) gradVLike = &m_tmpVector1;
103 M* hessianMLike = NULL;
104 if (hessianMatrix) hessianMLike = m_tmpMatrix;
106 V* hessianELike = NULL;
107 if (hessianEffect) hessianELike = &m_tmpVector2;
109 double value1 = m_priorDensity.actualValue (domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect);
111 if (m_likelihoodExponent != 0.) {
112 value2 = m_likelihoodFunction.actualValue(domainVector,domainDirection,gradVLike ,hessianMLike ,hessianELike );
117 "BayesianJointPdf<V,M>::actualValue()",
118 "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
120 double returnValue = value1;
121 if (m_likelihoodExponent == 0.) {
124 else if (m_likelihoodExponent == 1.) {
125 returnValue *= value2;
128 returnValue *= pow(value2,m_likelihoodExponent);
130 returnValue *= exp(m_logOfNormalizationFactor);
132 m_lastComputedLogPrior = log(value1);
133 m_lastComputedLogLikelihood = m_likelihoodExponent*log(value2);
135 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
136 *m_env.subDisplayFile() <<
"Leaving BayesianJointPdf<V,M>::actualValue()"
137 <<
": domainVector = " << domainVector
138 <<
", returnValue = " << returnValue
145 template<
class V,
class M>
148 const V& domainVector,
149 const V* domainDirection,
152 V* hessianEffect)
const
154 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
155 *m_env.subDisplayFile() <<
"Entering BayesianJointPdf<V,M>::lnValue()"
156 <<
": domainVector = " << domainVector
162 gradVLike = &m_tmpVector1;
167 for (
unsigned int i = 0; i < m_tmpVector1.sizeLocal(); i++) {
168 m_tmpVector1[i] = (*gradVector)[i];
172 M* hessianMLike = NULL;
173 if (hessianMatrix) hessianMLike = m_tmpMatrix;
175 V* hessianELike = NULL;
176 if (hessianEffect) hessianELike = &m_tmpVector2;
178 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
179 *m_env.subDisplayFile() <<
"In BayesianJointPdf<V,M>::lnValue()"
180 <<
", domainVector = " << domainVector
181 <<
": about to call prior()"
185 double value1 = m_priorDensity.lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect);
187 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
188 *m_env.subDisplayFile() <<
"In BayesianJointPdf<V,M>::lnValue()"
189 <<
", domainVector = " << domainVector
190 <<
": lnPrior = " << value1
194 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
195 *m_env.subDisplayFile() <<
"In BayesianJointPdf<V,M>::lnValue()"
196 <<
", domainVector = " << domainVector
197 <<
": about to call likelihood()"
202 if (m_likelihoodExponent != 0.) {
203 value2 = m_likelihoodFunction.lnValue(domainVector,domainDirection,gradVLike, hessianMLike, hessianELike );
206 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
207 *m_env.subDisplayFile() <<
"In BayesianJointPdf<V,M>::lnValue()"
208 <<
", domainVector = " << domainVector
209 <<
": value1 = " << value1
210 <<
", value2 = " << value2
213 *m_env.subDisplayFile() <<
"In BayesianJointPdf<V,M>::lnValue()"
214 <<
", domainVector = " << domainVector
215 <<
": gradVector = " << *gradVector
216 <<
", gradVLike = " << *gradVLike
220 *m_env.subDisplayFile() <<
"In BayesianJointPdf<V,M>::lnValue()"
221 <<
", domainVector = " << domainVector
222 <<
": hessianMatrix = " << *hessianMatrix
223 <<
", hessianMLike = " << *hessianMLike
227 *m_env.subDisplayFile() <<
"In BayesianJointPdf<V,M>::lnValue()"
228 <<
", domainVector = " << domainVector
229 <<
": hessianEffect = " << *hessianEffect
230 <<
", hessianELike = " << *hessianELike
235 if (gradVector ) *gradVector += *gradVLike;
236 if (hessianMatrix) *hessianMatrix += *hessianMLike;
237 if (hessianEffect) *hessianEffect += *hessianELike;
239 double returnValue = value1;
240 if (m_likelihoodExponent == 0.) {
243 else if (m_likelihoodExponent == 1.) {
244 returnValue += value2;
247 returnValue += value2*m_likelihoodExponent;
249 returnValue += m_logOfNormalizationFactor;
251 m_lastComputedLogPrior = value1;
252 m_lastComputedLogLikelihood = m_likelihoodExponent*value2;
254 if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
255 *m_env.subDisplayFile() <<
"Leaving BayesianJointPdf<V,M>::lnValue()"
256 <<
": domainVector = " << domainVector
257 <<
", returnValue = " << returnValue
264 template<
class V,
class M>
270 double volume = m_domainSet.volume();
271 if (((boost::math::isnan)(volume)) ||
272 (volume == -INFINITY ) ||
273 (volume == INFINITY ) ||
280 "BayesianJointPdf<V,M>::lnValue()",
281 "incomplete code for computeLogOfNormalizationFactor()");
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the PDF (scalar function).
double lastComputedLogLikelihood() const
Returns the logarithm of the last computed likelihood value. Access to protected attribute m_lastComp...
A templated class for handling sets.
~BayesianJointPdf()
Destructor.
double lnValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Computes the logarithm of the value of the function.
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
TODO: Computes the logarithm of the normalization factor.
void setNormalizationStyle(unsigned int value) const
Sets a value to be used in the normalization style of the prior density PDF (ie, protected attribute ...
A templated (base) class for handling scalar functions.
double lastComputedLogPrior() const
Returns the logarithm of the last computed Prior value. Access to protected attribute m_lastComputedL...
#define UQ_FATAL_TEST_MACRO(test, givenRank, where, what)
BayesianJointPdf(const char *prefix, const BaseJointPdf< V, M > &priorDensity, const BaseScalarFunction< V, M > &likelihoodFunction, double likelihoodExponent, const VectorSet< V, M > &intersectionDomain)
Default constructor.
A templated (base) class for handling joint PDFs.
A class for handling Bayesian joint PDFs.