queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
InverseGammaJointPdf.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/InverseGammaJointPdf.h>
26 #include <queso/GslVector.h>
27 #include <queso/GslMatrix.h>
28 
29 namespace QUESO {
30 
31 // Constructor -------------------------------------
32 template<class V,class M>
34  const char* prefix,
35  const VectorSet<V,M>& domainSet,
36  const V& alpha,
37  const V& beta)
38  :
39  BaseJointPdf<V,M>(((std::string)(prefix)+"uni").c_str(),domainSet),
40  m_alpha(alpha),
41  m_beta (beta)
42 {
43  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
44  *m_env.subDisplayFile() << "Entering InverseGammaJointPdf<V,M>::constructor()"
45  << ": prefix = " << m_prefix
46  << std::endl;
47  }
48 
49  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 54)) {
50  *m_env.subDisplayFile() << "Leaving InverseGammaJointPdf<V,M>::constructor()"
51  << ": prefix = " << m_prefix
52  << std::endl;
53  }
54 }
55 // Destructor --------------------------------------
56 template<class V,class M>
58 {
59 }
60 // Math methods-------------------------------------
61 template<class V, class M>
62 double
64  const V& domainVector,
65  const V* domainDirection,
66  V* gradVector,
67  M* hessianMatrix,
68  V* hessianEffect) const
69 {
70  queso_require_equal_to_msg(domainVector.sizeLocal(), this->m_domainSet.vectorSpace().dimLocal(), "invalid input");
71 
72  queso_require_msg(!(domainDirection || gradVector || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
73 
74  // No need to multiply by exp(m_logOfNormalizationFactor) because 'lnValue()' is called [PDF-07]
75  return exp(this->lnValue(domainVector,domainDirection,gradVector,hessianMatrix,hessianEffect));
76 }
77 //--------------------------------------------------
78 template<class V, class M>
79 double
81  const V& domainVector,
82  const V* domainDirection,
83  V* gradVector,
84  M* hessianMatrix,
85  V* hessianEffect) const
86 {
87  queso_require_msg(!(domainDirection || gradVector || hessianMatrix || hessianEffect), "incomplete code for gradVector, hessianMatrix and hessianEffect calculations");
88 
89  double result = 0.;
90  for (unsigned int i = 0; i < domainVector.sizeLocal(); ++i) {
91  result -= (m_alpha[i]+1.)*log(domainVector[i]);
92  result -= m_beta[i]/domainVector[i];
93  if (m_normalizationStyle == 0) {
94  // Code needs to be done yet
95  }
96  }
97  result += m_logOfNormalizationFactor; // [PDF-07]
98 
99  return result;
100 }
101 //--------------------------------------------------
102 template<class V, class M>
103 void
105 {
106  queso_assert_equal_to(m_alpha.sizeLocal(), m_beta.sizeLocal());
107  queso_assert_equal_to(m_alpha.sizeLocal(), meanVector.sizeLocal());
108 
109  for (unsigned int i = 0; i < m_alpha.sizeLocal(); ++i) {
110  queso_assert_greater(m_alpha[i], 1);
111  meanVector[i] = m_beta[i] / (m_alpha[i] - 1);
112  }
113 }
114 //--------------------------------------------------
115 template<class V, class M>
116 void
118 {
119  queso_assert_equal_to(m_alpha.sizeLocal(), m_beta.sizeLocal());
120  queso_assert_equal_to(m_alpha.sizeLocal(), covMatrix.numCols());
121  queso_assert_equal_to (covMatrix.numCols(), covMatrix.numRowsGlobal());
122 
123  covMatrix.zeroLower();
124  covMatrix.zeroUpper();
125 
126  for (unsigned int i = 0; i < m_alpha.sizeLocal(); ++i) {
127  queso_assert_greater(m_alpha[i], 2);
128  covMatrix(i,i) = m_beta[i]*m_beta[i] / (m_alpha[i] - 1) /
129  (m_alpha[i] - 1) / (m_alpha[i] - 2);
130  }
131 }
132 //--------------------------------------------------
133 template<class V, class M>
134 double
135 InverseGammaJointPdf<V,M>::computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
136 {
137  double value = 0.;
138 
139  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
140  *m_env.subDisplayFile() << "Entering InverseGammaJointPdf<V,M>::computeLogOfNormalizationFactor()"
141  << std::endl;
142  }
143  value = BaseJointPdf<V,M>::commonComputeLogOfNormalizationFactor(numSamples, updateFactorInternally);
144  if ((m_env.subDisplayFile()) && (m_env.displayVerbosity() >= 2)) {
145  *m_env.subDisplayFile() << "Leaving InverseGammaJointPdf<V,M>::computeLogOfNormalizationFactor()"
146  << ", m_logOfNormalizationFactor = " << m_logOfNormalizationFactor
147  << std::endl;
148  }
149 
150  return value;
151 }
152 
153 } // End namespace QUESO
154 
double actualValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
Actual value of the Gamma PDF.
double lnValue(const V &domainVector, const V *domainDirection, V *gradVector, M *hessianMatrix, V *hessianEffect) const
TODO: Logarithm of the value of the Gamma PDF.
virtual void distributionVariance(M &covMatrix) const
Covariance matrix of the underlying random variable.
A templated class for handling sets.
Definition: VectorSet.h:52
double computeLogOfNormalizationFactor(unsigned int numSamples, bool updateFactorInternally) const
Computes the logarithm of the normalization factor.
InverseGammaJointPdf(const char *prefix, const VectorSet< V, M > &domainSet, const V &alpha, const V &beta)
Constructor.
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 class for handling Inverse Gamma joint PDFs.
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"))
unsigned int displayVerbosity() const
Definition: Environment.C:450
virtual void distributionMean(V &meanVector) const
Mean value of the underlying random variable.
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