queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
BetaJointPdf.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #include <queso/BetaJointPdf.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
28 #include <queso/BasicPdfsBase.h>
29 
30 namespace QUESO {
31 
32 // Constructor -------------------------------------
33 template<class V,class M>
35  const char* prefix,
36  const VectorSet<V,M>& domainSet,
37  const V& alpha,
38  const V& beta)
39  :
40  BaseJointPdf<V,M>(((std::string)(prefix)+"uni").c_str(),domainSet),
41  m_alpha(alpha),
42  m_beta (beta)
43 {
44  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
45  *m_env.subDisplayFile() << "Entering BetaJointPdf<V,M>::constructor()"
46  << ": prefix = " << m_prefix
47  << std::endl;
48  }
49 
50  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
51  *m_env.subDisplayFile() << "Leaving BetaJointPdf<V,M>::constructor()"
52  << ": prefix = " << m_prefix
53  << std::endl;
54  }
55 }
56 // Destructor --------------------------------------
57 template<class V,class M>
59 {
60 }
61 // Math methods-------------------------------------
62 template<class V, class M>
63 double
65  const V& domainVector,
66  const V* domainDirection,
67  V* gradVector,
68  M* hessianMatrix,
69  V* hessianEffect) const
70 {
71  queso_require_equal_to_msg(domainVector.sizeLocal(), this->m_domainSet.vectorSpace().dimLocal(), "invalid input");
72  queso_require_msg(!(domainDirection || hessianMatrix || hessianEffect), "incomplete code for hessianMatrix and hessianEffect calculations");
73 
74  // No need to multiply by exp(m_logOfNormalizationFactor) because 'lnValue()' is called
75  double value = std::exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
76 
77  if (gradVector) {
78  // DM: This transformation may have issues near the boundary. I haven't
79  // fully explored this yet.
80  //
81  // This evaluation is also normalised (if m_normalizationStyle is zero) and
82  // so the gradient in physical space contains the normalisation constant.
83  (*gradVector) *= value;
84  }
85 
86  return value;
87 }
88 //--------------------------------------------------
89 template<class V, class M>
90 double
92  const V& domainVector,
93  const V* domainDirection,
94  V* gradVector,
95  M* hessianMatrix,
96  V* hessianEffect) const
97 {
98  queso_require_msg(!(domainDirection || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
99 
100  double aux = 0.;
101  double result = 0.;
102  for (unsigned int i = 0; i < domainVector.sizeLocal(); ++i) {
103  if (m_normalizationStyle == 0) {
104  aux = log(m_env.basicPdfs()->betaPdfActualValue(domainVector[i],m_alpha[i],m_beta[i]));
105  }
106  else {
107  aux = (m_alpha[i]-1.)*log(domainVector[i]) + (m_beta[i]-1.)*log(1.-domainVector[i]);
108  }
109 
110  // The log of the beta pdf is
111  // f(x) = (alpha - 1) * log(x) + (beta - 1) * log(1 - x)
112  // Therefore
113  // df/dx (x) = ((alpha - 1) / x) + ((1 - beta) / (1 - x))
114  if (gradVector) {
115  // We're computing grad of log of p which is p' / p
116  (*gradVector)[i] = (m_alpha[i] - 1.0) / domainVector[i] +
117  (1.0 - m_beta[i]) / (1.0 - domainVector[i]);
118  }
119 
120  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 99)) {
121  *m_env.subDisplayFile() << "In BetaJointPdf<V,M>::lnValue()"
122  << ", m_normalizationStyle = " << m_normalizationStyle
123  << ": domainVector[" << i << "] = " << domainVector[i]
124  << ", m_alpha[" << i << "] = " << m_alpha[i]
125  << ", m_beta[" << i << "] = " << m_beta[i]
126  << ", log(pdf)= " << aux
127  << std::endl;
128  }
129  result += aux;
130  }
131  result += m_logOfNormalizationFactor;
132 
133  return result;
134 }
135 //--------------------------------------------------
136 template<class V, class M>
137 void
139 {
140  unsigned int n_params = meanVector.sizeLocal();
141  queso_assert_equal_to (n_params, m_alpha.sizeLocal());
142 
143  for (unsigned int i = 0; i < n_params; ++i) {
144  meanVector[i] = m_alpha[i] / (m_alpha[i] + m_beta[i]);
145  }
146 }
147 //--------------------------------------------------
148 template<class V, class M>
149 void
151 {
152  unsigned int n_params = m_alpha.sizeLocal();
153  queso_assert_equal_to (n_params, m_beta.sizeLocal());
154  queso_assert_equal_to (n_params, covMatrix.numCols());
155  queso_assert_equal_to (covMatrix.numCols(), covMatrix.numRowsGlobal());
156 
157  covMatrix.zeroLower();
158  covMatrix.zeroUpper();
159 
160  for (unsigned int i = 0; i < n_params; ++i) {
161  covMatrix(i,i) = (m_alpha[i] * m_beta[i]) /
162  ((m_alpha[i] + m_beta[i]) *
163  (m_alpha[i] + m_beta[i]) *
164  (m_alpha[i] + m_beta[i] + 1));
165  }
166 }
167 //--------------------------------------------------
168 template<class V, class M>
169 double
170 BetaJointPdf<V,M>::computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
171 {
172  double value = 0.;
173 
174  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
175  *m_env.subDisplayFile() << "Entering BetaJointPdf<V,M>::computeLogOfNormalizationFactor()"
176  << std::endl;
177  }
178  value = BaseJointPdf<V,M>::commonComputeLogOfNormalizationFactor(numSamples, updateFactorInternally);
179  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
180  *m_env.subDisplayFile() << "Leaving BetaJointPdf<V,M>::computeLogOfNormalizationFactor()"
181  << ", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
182  << std::endl;
183  }
184 
185  return value;
186 }
187 
188 } // End namespace QUESO
189 
BetaJointPdf(const char *prefix, const VectorSet< V, M > &domainSet, const V &alpha, const V &beta)
Constructor.
Definition: BetaJointPdf.C:34
double lnValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Logarithm of the value of the Beta PDF.
Definition: BetaJointPdf.C:91
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
Definition: BetaJointPdf.C:170
virtual void distributionVariance(M &covMatrix) const
Covariance matrix of the underlying random variable.
Definition: BetaJointPdf.C:150
A templated class for handling sets.
Definition: VectorSet.h:52
~BetaJointPdf()
Destructor.
Definition: BetaJointPdf.C:58
double commonComputeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Common method (to the derived classes) to compute the logarithm of the normalization factor...
Definition: JointPdf.C:115
const BaseEnvironment & m_env
A templated (base) class for handling joint PDFs.
Definition: JointPdf.h:55
MonteCarloSGOptions::MonteCarloSGOptions(const BaseEnvironment &env, const char *prefix, const McOptionsValues &alternativeOptionsValues queso_require_equal_to_msg)(m_env.optionsInputFileName(), std::string(""), std::string("this constructor is incompatible with the existence of an options input file"))
virtual void distributionMean(V &meanVector) const
Mean value of the underlying random variable.
Definition: BetaJointPdf.C:138
unsigned int displayVerbosity() const
Definition: Environment.C:450
A class for handling Beta joint PDFs.
Definition: BetaJointPdf.h:47
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the Beta PDF.
Definition: BetaJointPdf.C:64
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320

Generated on Tue Jun 5 2018 19:48:55 for queso-0.57.1 by  doxygen 1.8.5