From 555612d2a6f750321d2e76d715b49cab5051f6e4 Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Mon, 6 Feb 2017 11:10:40 -0500 Subject: [PATCH 1/8] Remove old debugging print statement from test --- test/unit/quadrature_1d.C | 1 - 1 file changed, 1 deletion(-) diff --git a/test/unit/quadrature_1d.C b/test/unit/quadrature_1d.C index 75265ebe6..82826d7b5 100644 --- a/test/unit/quadrature_1d.C +++ b/test/unit/quadrature_1d.C @@ -166,7 +166,6 @@ namespace QUESOTesting void test_1d_hermite_quadrature_x4erf() { - std::cout << std::endl << "Beginning Hermite X4Erf test!" << std::endl; X4Erf func; std::vector testing_orders; this->testing_orders(testing_orders); From 3136a4f5acc52a609fcff6d811a332788d77b27e Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Tue, 24 Jan 2017 13:18:32 -0500 Subject: [PATCH 2/8] Introduce LikehoodBase::lnLikelihood() We want to factor out the model evaluation from the actual computation of the log-likelihood for forthcoming additions of "marginalized parameters". So, we add lnValue() to LikelihoodBase, where it calls evaluateModel() first, then we pass the model evaluation to lnLikelihood(). Now, subclasses implement lnLikelihood with the modelOutput as an input parameter to the function. We don't pass modelOutput as const so that it can be overwritten to avoid copying the modelOutput, preserving the behavior that was there before. I'm not sure I like this behavior, but it's also a protected function that's working a vector created internal to the class hierarchy. --- .../inc/GaussianLikelihoodBlockDiagonalCovariance.h | 6 ++++-- ...lihoodBlockDiagonalCovarianceRandomCoefficients.h | 10 ++++------ src/stats/inc/GaussianLikelihoodDiagonalCovariance.h | 6 ++++-- src/stats/inc/GaussianLikelihoodFullCovariance.h | 6 ++++-- ...ussianLikelihoodFullCovarianceRandomCoefficient.h | 6 ++++-- src/stats/inc/GaussianLikelihoodScalarCovariance.h | 6 ++++-- src/stats/inc/LikelihoodBase.h | 12 ++++++++++++ .../src/GaussianLikelihoodBlockDiagonalCovariance.C | 5 +---- ...lihoodBlockDiagonalCovarianceRandomCoefficients.C | 6 ++---- src/stats/src/GaussianLikelihoodDiagonalCovariance.C | 6 +----- src/stats/src/GaussianLikelihoodFullCovariance.C | 5 +---- ...ussianLikelihoodFullCovarianceRandomCoefficient.C | 5 +---- src/stats/src/GaussianLikelihoodScalarCovariance.C | 6 +----- src/stats/src/LikelihoodBase.C | 11 +++++++++++ 14 files changed, 54 insertions(+), 42 deletions(-) diff --git a/src/stats/inc/GaussianLikelihoodBlockDiagonalCovariance.h b/src/stats/inc/GaussianLikelihoodBlockDiagonalCovariance.h index cc0ae5422..44db646bb 100644 --- a/src/stats/inc/GaussianLikelihoodBlockDiagonalCovariance.h +++ b/src/stats/inc/GaussianLikelihoodBlockDiagonalCovariance.h @@ -70,8 +70,10 @@ class GaussianLikelihoodBlockDiagonalCovariance : public LikelihoodBase { //! Get (const) multiplicative coefficient for block \c i const double & getBlockCoefficient(unsigned int i) const; - //! Logarithm of the value of the scalar function. - virtual double lnValue(const V & domainVector) const; +protected: + + //! Logarithm of the value of the likelihood function. + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const; private: std::vector m_covarianceCoefficients; diff --git a/src/stats/inc/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h b/src/stats/inc/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h index afee9e955..74f6b55cd 100644 --- a/src/stats/inc/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h +++ b/src/stats/inc/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h @@ -66,17 +66,15 @@ class GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients : public Likel virtual ~GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients(); //@} - //! Logarithm of the value of the scalar function. +protected: + + //! Logarithm of the value of the likelihood function. /*! * The last \c n elements of \c domainVector, where \c n is the number of * blocks in the block diagonal covariance matrix, are treated as * hyperparameters and will be sample in a statistical inversion. - * - * The user need not concern themselves with handling these in the model - * evaluation routine, since they are handled by the likelihood evaluation - * routine. */ - virtual double lnValue(const V & domainVector) const; + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const; private: const GslBlockMatrix & m_covariance; diff --git a/src/stats/inc/GaussianLikelihoodDiagonalCovariance.h b/src/stats/inc/GaussianLikelihoodDiagonalCovariance.h index 4c291a079..438bf2601 100644 --- a/src/stats/inc/GaussianLikelihoodDiagonalCovariance.h +++ b/src/stats/inc/GaussianLikelihoodDiagonalCovariance.h @@ -58,8 +58,10 @@ class GaussianLikelihoodDiagonalCovariance : public LikelihoodBase { virtual ~GaussianLikelihoodDiagonalCovariance(); //@} - //! Logarithm of the value of the scalar function. - virtual double lnValue(const V & domainVector) const; +protected: + + //! Logarithm of the value of the likelihood function. + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const; private: const V & m_covariance; diff --git a/src/stats/inc/GaussianLikelihoodFullCovariance.h b/src/stats/inc/GaussianLikelihoodFullCovariance.h index ff6edbf43..1986d505a 100644 --- a/src/stats/inc/GaussianLikelihoodFullCovariance.h +++ b/src/stats/inc/GaussianLikelihoodFullCovariance.h @@ -62,8 +62,10 @@ class GaussianLikelihoodFullCovariance : public LikelihoodBase { virtual ~GaussianLikelihoodFullCovariance(); //@} - //! Logarithm of the value of the scalar function. - virtual double lnValue(const V & domainVector) const; +protected: + + //! Logarithm of the value of the likelihood function. + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const; private: double m_covarianceCoefficient; diff --git a/src/stats/inc/GaussianLikelihoodFullCovarianceRandomCoefficient.h b/src/stats/inc/GaussianLikelihoodFullCovarianceRandomCoefficient.h index 2ed5c7ef1..9a71b96f6 100644 --- a/src/stats/inc/GaussianLikelihoodFullCovarianceRandomCoefficient.h +++ b/src/stats/inc/GaussianLikelihoodFullCovarianceRandomCoefficient.h @@ -67,8 +67,10 @@ class GaussianLikelihoodFullCovarianceRandomCoefficient : public LikelihoodBase< virtual ~GaussianLikelihoodFullCovarianceRandomCoefficient(); //@} - //! Logarithm of the value of the scalar function. - virtual double lnValue(const V & domainVector) const; +protected: + + //! Logarithm of the value of the likelihood function. + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const; private: const M & m_covariance; diff --git a/src/stats/inc/GaussianLikelihoodScalarCovariance.h b/src/stats/inc/GaussianLikelihoodScalarCovariance.h index ab5a459bd..ec1c072e1 100644 --- a/src/stats/inc/GaussianLikelihoodScalarCovariance.h +++ b/src/stats/inc/GaussianLikelihoodScalarCovariance.h @@ -58,8 +58,10 @@ class GaussianLikelihoodScalarCovariance : public LikelihoodBase { virtual ~GaussianLikelihoodScalarCovariance(); //@} - //! Logarithm of the value of the scalar function. - virtual double lnValue(const V & domainVector) const; +protected: + + //! Logarithm of the value of the likelihood function. + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const; private: double m_covariance; diff --git a/src/stats/inc/LikelihoodBase.h b/src/stats/inc/LikelihoodBase.h index fb53754e3..68567d8cc 100644 --- a/src/stats/inc/LikelihoodBase.h +++ b/src/stats/inc/LikelihoodBase.h @@ -99,8 +99,20 @@ class LikelihoodBase : public BaseScalarFunction { V * /*gradVector*/, M * /*hessianMatrix*/, V * /*hessianEffect*/) const { return std::exp(this->lnValue(domainVector)); } + virtual double lnValue(const V & domainVector) const; + protected: + const V & m_observations; + + //! Compute log-likelihood value given the current parameters and model output + /*! + * Subclasses should override this method to compute the likelihood value, given + * the current parameters in domainVector and the output from the model evaluated + * at those parameter values. Note that this function may alter the values of + * modelOutput to reduce copying. + */ + virtual double lnLikelihood(const V & domainVector, V & modelOutput) const =0; }; } // End namespace QUESO diff --git a/src/stats/src/GaussianLikelihoodBlockDiagonalCovariance.C b/src/stats/src/GaussianLikelihoodBlockDiagonalCovariance.C index ee3813a30..9b8f9f62c 100644 --- a/src/stats/src/GaussianLikelihoodBlockDiagonalCovariance.C +++ b/src/stats/src/GaussianLikelihoodBlockDiagonalCovariance.C @@ -73,13 +73,10 @@ GaussianLikelihoodBlockDiagonalCovariance::getBlockCoefficient( template double -GaussianLikelihoodBlockDiagonalCovariance::lnValue(const V & domainVector) const +GaussianLikelihoodBlockDiagonalCovariance::lnLikelihood(const V & domainVector, V & modelOutput) const { - V modelOutput(this->m_observations, 0, 0); // At least it's not a copy V weightedMisfit(this->m_observations, 0, 0); // At least it's not a copy - this->evaluateModel(domainVector, modelOutput); - // Compute misfit G(x) - y modelOutput -= this->m_observations; diff --git a/src/stats/src/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.C b/src/stats/src/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.C index 3933392fb..a84c6ca0f 100644 --- a/src/stats/src/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.C +++ b/src/stats/src/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.C @@ -54,13 +54,11 @@ GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::~GaussianLike template double -GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::lnValue(const V & domainVector) const +GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::lnLikelihood(const V & domainVector, + V & modelOutput) const { - V modelOutput(this->m_observations, 0, 0); // At least it's not a copy V weightedMisfit(this->m_observations, 0, 0); // At least it's not a copy - this->evaluateModel(domainVector, modelOutput); - // Compute misfit G(x) - y modelOutput -= this->m_observations; diff --git a/src/stats/src/GaussianLikelihoodDiagonalCovariance.C b/src/stats/src/GaussianLikelihoodDiagonalCovariance.C index 639e45d1d..0c94bd75f 100644 --- a/src/stats/src/GaussianLikelihoodDiagonalCovariance.C +++ b/src/stats/src/GaussianLikelihoodDiagonalCovariance.C @@ -50,12 +50,8 @@ GaussianLikelihoodDiagonalCovariance::~GaussianLikelihoodDiagonalCovarianc template double -GaussianLikelihoodDiagonalCovariance::lnValue(const V & domainVector) const +GaussianLikelihoodDiagonalCovariance::lnLikelihood(const V & domainVector, V & modelOutput) const { - V modelOutput(this->m_observations, 0, 0); // At least it's not a copy - - this->evaluateModel(domainVector, modelOutput); - modelOutput -= this->m_observations; // Compute misfit modelOutput *= modelOutput; modelOutput /= this->m_covariance; // Multiply by inverse covriance matrix diff --git a/src/stats/src/GaussianLikelihoodFullCovariance.C b/src/stats/src/GaussianLikelihoodFullCovariance.C index 3e8a0bd2f..090f4e2b3 100644 --- a/src/stats/src/GaussianLikelihoodFullCovariance.C +++ b/src/stats/src/GaussianLikelihoodFullCovariance.C @@ -51,13 +51,10 @@ GaussianLikelihoodFullCovariance::~GaussianLikelihoodFullCovariance() template double -GaussianLikelihoodFullCovariance::lnValue(const V & domainVector) const +GaussianLikelihoodFullCovariance::lnLikelihood(const V & domainVector, V & modelOutput) const { - V modelOutput(this->m_observations, 0, 0); // At least it's not a copy V weightedMisfit(this->m_observations, 0, 0); // At least it's not a copy - this->evaluateModel(domainVector, modelOutput); - // Compute misfit G(x) - y modelOutput -= this->m_observations; diff --git a/src/stats/src/GaussianLikelihoodFullCovarianceRandomCoefficient.C b/src/stats/src/GaussianLikelihoodFullCovarianceRandomCoefficient.C index 349c51a5c..65481b12e 100644 --- a/src/stats/src/GaussianLikelihoodFullCovarianceRandomCoefficient.C +++ b/src/stats/src/GaussianLikelihoodFullCovarianceRandomCoefficient.C @@ -49,13 +49,10 @@ GaussianLikelihoodFullCovarianceRandomCoefficient::~GaussianLikelihoodFull template double -GaussianLikelihoodFullCovarianceRandomCoefficient::lnValue(const V & domainVector) const +GaussianLikelihoodFullCovarianceRandomCoefficient::lnLikelihood(const V & domainVector, V & modelOutput) const { - V modelOutput(this->m_observations, 0, 0); // At least it's not a copy V weightedMisfit(this->m_observations, 0, 0); // At least it's not a copy - this->evaluateModel(domainVector, modelOutput); - // Compute misfit G(x) - y modelOutput -= this->m_observations; diff --git a/src/stats/src/GaussianLikelihoodScalarCovariance.C b/src/stats/src/GaussianLikelihoodScalarCovariance.C index a6ca88e82..c74ab1f93 100644 --- a/src/stats/src/GaussianLikelihoodScalarCovariance.C +++ b/src/stats/src/GaussianLikelihoodScalarCovariance.C @@ -47,12 +47,8 @@ GaussianLikelihoodScalarCovariance::~GaussianLikelihoodScalarCovariance() template double -GaussianLikelihoodScalarCovariance::lnValue(const V & domainVector) const +GaussianLikelihoodScalarCovariance::lnLikelihood(const V & domainVector, V & modelOutput) const { - V modelOutput(this->m_observations, 0, 0); // At least it's not a copy - - this->evaluateModel(domainVector, modelOutput); - modelOutput -= this->m_observations; // Compute misfit double norm2_squared = modelOutput.norm2Sq(); // Compute square of 2-norm diff --git a/src/stats/src/LikelihoodBase.C b/src/stats/src/LikelihoodBase.C index f81a0271e..1fc4d342b 100644 --- a/src/stats/src/LikelihoodBase.C +++ b/src/stats/src/LikelihoodBase.C @@ -59,6 +59,17 @@ void LikelihoodBase::evaluateModel(const V & domainVector, const V * domai queso_error_msg(ss.str()); } +template +double +LikelihoodBase::lnValue(const V & domainVector) const +{ + V modelOutput(this->m_observations, 0, 0); // At least it's not a copy + + this->evaluateModel(domainVector, modelOutput); + + return this->lnLikelihood(domainVector, modelOutput); +} + } // End namespace QUESO template class QUESO::LikelihoodBase; From 83717003bcb41af7a7bd37bb4168ca7cc1c467f5 Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Tue, 24 Jan 2017 17:09:47 -0500 Subject: [PATCH 3/8] Add new LikelihoodBase constructor for marginal params This is laying the groundwork for having the likelihood do marginalization for the user. This add the neccessary data and constructor to LikelihoodBase. In particular, we'll need the joint pdf of the marginal parameter space and the user-constructed integration rule for integrating over that space. --- src/stats/inc/LikelihoodBase.h | 45 ++++++++++++++++++++++++++++++++++ src/stats/src/LikelihoodBase.C | 21 ++++++++++++++++ 2 files changed, 66 insertions(+) diff --git a/src/stats/inc/LikelihoodBase.h b/src/stats/inc/LikelihoodBase.h index 68567d8cc..991c55b73 100644 --- a/src/stats/inc/LikelihoodBase.h +++ b/src/stats/inc/LikelihoodBase.h @@ -28,9 +28,16 @@ #include #include #include +#include namespace QUESO { +template +class BaseVectorRV; + +template +class MultiDQuadratureBase; + class GslVector; class GslMatrix; @@ -59,8 +66,34 @@ class LikelihoodBase : public BaseScalarFunction { const VectorSet & domainSet, const V & observations); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + LikelihoodBase(const char * prefix, + const VectorSet & domainSet, + const V & observations, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor, pure to make this class abstract virtual ~LikelihoodBase() =0; + //@} //! Deprecated. Evaluates the user's model at the point \c domainVector @@ -99,12 +132,24 @@ class LikelihoodBase : public BaseScalarFunction { V * /*gradVector*/, M * /*hessianMatrix*/, V * /*hessianEffect*/) const { return std::exp(this->lnValue(domainVector)); } + //! Logarithm of the value of the scalar function. + /*! + * This method well evaluate the users model (evaluateModel()) and then + * call lnLikelihood() and, if there is marginalization, will handle the + * numerical integration over the marginal parameter space. + */ virtual double lnValue(const V & domainVector) const; protected: const V & m_observations; + typename SharedPtr >::Type m_marg_param_pdf; + + typename SharedPtr >::Type m_marg_integration; + + bool m_marg_pdf_is_weight_func; + //! Compute log-likelihood value given the current parameters and model output /*! * Subclasses should override this method to compute the likelihood value, given diff --git a/src/stats/src/LikelihoodBase.C b/src/stats/src/LikelihoodBase.C index 1fc4d342b..d896f4c3a 100644 --- a/src/stats/src/LikelihoodBase.C +++ b/src/stats/src/LikelihoodBase.C @@ -28,6 +28,8 @@ #include #include #include +#include +#include namespace QUESO { @@ -40,6 +42,25 @@ LikelihoodBase::LikelihoodBase( { } +template +LikelihoodBase::LikelihoodBase(const char * prefix, + const VectorSet & domainSet, + const V & observations, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : BaseScalarFunction(prefix, domainSet), + m_observations(observations), + m_marg_param_pdf(marg_param_pdf), + m_marg_integration(marg_integration), + m_marg_pdf_is_weight_func(marg_pdf_is_weight_func) +{ + // The dimension of the parameter space had better match the dimension in the integration + queso_require_equal_to_msg(this->m_marg_param_pdf->imageSet().vectorSpace().dimGlobal(), + this->m_marg_integration->getDomain().vectorSpace().dimGlobal(), + "Mismatched marginal parameter space dimension and quadrature dimension!"); +} + template LikelihoodBase::~LikelihoodBase() { From 0d077b48b3e7bb66011eedd009e5c195651ed80d Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Tue, 24 Jan 2017 17:19:31 -0500 Subject: [PATCH 4/8] Add new evaluateModel API for likelihood marginalization The user will need to override this method for doing marginalization in the likelihood. In particular, this method will also pass the current value of the marginal parameters in addition to the inference parameters. This will be the user function that gets called if a marginal parameter space and integration rule have been provided. --- src/stats/inc/LikelihoodBase.h | 17 ++++++++++++++++- src/stats/src/LikelihoodBase.C | 14 ++++++++++++++ 2 files changed, 30 insertions(+), 1 deletion(-) diff --git a/src/stats/inc/LikelihoodBase.h b/src/stats/inc/LikelihoodBase.h index 991c55b73..20a39286c 100644 --- a/src/stats/inc/LikelihoodBase.h +++ b/src/stats/inc/LikelihoodBase.h @@ -112,7 +112,8 @@ class LikelihoodBase : public BaseScalarFunction { * be compared to actual observations when computing the likelihood functional. * * The first \c n components of \c domainVector are the model parameters. - * The rest of \c domainVector contains the hyperparameters, if any. For + * The rest of \c domainVector contains the hyperparameters, if any. Mathematically, this + * function is \f$ f(m) \f$. For * example, in \c GaussianLikelihoodFullCovarainceRandomCoefficient, the last * component of \c domainVector contains the multiplicative coefficient of * the observational covariance matrix. In this case, the user need not @@ -127,6 +128,20 @@ class LikelihoodBase : public BaseScalarFunction { virtual void evaluateModel(const V & domainVector, V & modelOutput) const { this->evaluateModel(domainVector,NULL,modelOutput,NULL,NULL,NULL); } + //! Evaluates the user's model at the point \c (domainVector,marginalVector) + /*! + * Subclass implementations will fill up the \c modelOutput vector with output + * from the model. This function will be passed the current value of the parameters in domainVector + * and the current values of the marginalization parameters in marginalVector. Mathematically, this + * function is \f$ f(m,q) \f$. This function + * is only called if the likelihood includes marginalization. Otherwise, + * evaluateModel(const V & domainVector, V & modelOutput) will be called by lnValue(). + * + * Not every user will do marginalization, so this is not a pure function, but there is + * no default, so we error if marginalization is attempted without overriding this function. + */ + virtual void evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput) const; + //! Actual value of the scalar function. virtual double actualValue(const V & domainVector, const V * /*domainDirection*/, V * /*gradVector*/, M * /*hessianMatrix*/, V * /*hessianEffect*/) const diff --git a/src/stats/src/LikelihoodBase.C b/src/stats/src/LikelihoodBase.C index d896f4c3a..6b09722ad 100644 --- a/src/stats/src/LikelihoodBase.C +++ b/src/stats/src/LikelihoodBase.C @@ -80,6 +80,20 @@ void LikelihoodBase::evaluateModel(const V & domainVector, const V * domai queso_error_msg(ss.str()); } +template +void LikelihoodBase::evaluateModel(const V & domainVector, + const V & marginalVector, + V & modelOutput) const +{ + std::stringstream ss; + ss << "ERROR: evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput)" + << std::endl + << " is not implemented! Please override this function in your subclass." + << std::endl; + + queso_error_msg(ss.str()); +} + template double LikelihoodBase::lnValue(const V & domainVector) const From a52f42d84d07055d0527a50b74ee8119f7740f1d Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Tue, 24 Jan 2017 17:21:45 -0500 Subject: [PATCH 5/8] Likelihood evaluation can now do marginalization If we detect the user has supplied a marinal parameter pdf, then we instead will use the user-supplied integration rule to integrate over the parameter space, calling the user-overrideen evaluateModel() method that accepts marginal parameters as well as inference parameters. --- src/stats/src/LikelihoodBase.C | 49 +++++++++++++++++++++++++++++++--- 1 file changed, 46 insertions(+), 3 deletions(-) diff --git a/src/stats/src/LikelihoodBase.C b/src/stats/src/LikelihoodBase.C index 6b09722ad..d45ff1a0e 100644 --- a/src/stats/src/LikelihoodBase.C +++ b/src/stats/src/LikelihoodBase.C @@ -98,11 +98,54 @@ template double LikelihoodBase::lnValue(const V & domainVector) const { - V modelOutput(this->m_observations, 0, 0); // At least it's not a copy + V modelOutput(this->m_observations, 0, 0); - this->evaluateModel(domainVector, modelOutput); + double lnLikelihood_value = 0.0; - return this->lnLikelihood(domainVector, modelOutput); + // If we're not marginalizing, then we just need f(m) + if( !m_marg_param_pdf ) + { + this->evaluateModel(domainVector, modelOutput); + + // Note modelOutput made be modified in lnLikelihood() + lnLikelihood_value = this->lnLikelihood(domainVector, modelOutput); + } + // Otherwise we're integrating over marginal parameter space + else + { + queso_assert(m_marg_integration); + + const std::vector< typename QUESO::SharedPtr< V >::Type > & x = + this->m_marg_integration->positions(); + + const std::vector< double > & w = this->m_marg_integration->weights(); + + unsigned int n_qpoints = x.size(); + + for( unsigned int q = 0; q < n_qpoints; q++ ) + { + this->evaluateModel(domainVector, *(x[q]), modelOutput); + + // Note modelOutput made be modified in lnLikelihood() + double ln_pi_m_q = this->lnLikelihood(domainVector, modelOutput); + + double ln_pi_q = 0.0; + + // If we're not treating the marginal parameters pdf as + // the weighting function in the quadrature evaluation, + // then we need to also evaluate it. + if(!m_marg_pdf_is_weight_func) + this->m_marg_param_pdf->pdf().lnValue( *(x[q]) ); + + /*! \todo [PB]: We might want to play games with a log(\sum) identity + if precision starts to become an issue. */ + lnLikelihood_value += std::exp(ln_pi_m_q + ln_pi_q)*w[q]; + } + + lnLikelihood_value = std::log(lnLikelihood_value); + } + + return lnLikelihood_value; } } // End namespace QUESO From 5c0d380fee52696b8a5fd5264e7fdd104ec95580 Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Fri, 10 Feb 2017 09:04:45 -0500 Subject: [PATCH 6/8] Add marginalization constructor to Gaussian subclasses --- ...aussianLikelihoodBlockDiagonalCovariance.h | 28 ++++++++++++++ ...lockDiagonalCovarianceRandomCoefficients.h | 28 ++++++++++++++ .../GaussianLikelihoodDiagonalCovariance.h | 26 +++++++++++++ .../inc/GaussianLikelihoodFullCovariance.h | 25 ++++++++++++ ...ikelihoodFullCovarianceRandomCoefficient.h | 25 ++++++++++++ .../inc/GaussianLikelihoodScalarCovariance.h | 25 ++++++++++++ ...aussianLikelihoodBlockDiagonalCovariance.C | 38 +++++++++++++++---- ...lockDiagonalCovarianceRandomCoefficients.C | 37 ++++++++++++++---- .../GaussianLikelihoodDiagonalCovariance.C | 15 ++++++++ .../src/GaussianLikelihoodFullCovariance.C | 16 ++++++++ ...ikelihoodFullCovarianceRandomCoefficient.C | 15 ++++++++ .../src/GaussianLikelihoodScalarCovariance.C | 12 ++++++ 12 files changed, 274 insertions(+), 16 deletions(-) diff --git a/src/stats/inc/GaussianLikelihoodBlockDiagonalCovariance.h b/src/stats/inc/GaussianLikelihoodBlockDiagonalCovariance.h index 44db646bb..7c2e46476 100644 --- a/src/stats/inc/GaussianLikelihoodBlockDiagonalCovariance.h +++ b/src/stats/inc/GaussianLikelihoodBlockDiagonalCovariance.h @@ -60,6 +60,32 @@ class GaussianLikelihoodBlockDiagonalCovariance : public LikelihoodBase { const VectorSet & domainSet, const V & observations, const GslBlockMatrix & covariance); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + GaussianLikelihoodBlockDiagonalCovariance(const char * prefix, + const VectorSet & domainSet, + const V & observations, + const GslBlockMatrix & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor virtual ~GaussianLikelihoodBlockDiagonalCovariance(); //@} @@ -78,6 +104,8 @@ class GaussianLikelihoodBlockDiagonalCovariance : public LikelihoodBase { private: std::vector m_covarianceCoefficients; const GslBlockMatrix & m_covariance; + + void checkDimConsistency(const V & observations, const GslBlockMatrix & covariance) const; }; } // End namespace QUESO diff --git a/src/stats/inc/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h b/src/stats/inc/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h index 74f6b55cd..426863092 100644 --- a/src/stats/inc/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h +++ b/src/stats/inc/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h @@ -62,6 +62,32 @@ class GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients : public Likel const VectorSet & domainSet, const V & observations, const GslBlockMatrix & covariance); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients(const char * prefix, + const VectorSet & domainSet, + const V & observations, + const GslBlockMatrix & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor virtual ~GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients(); //@} @@ -78,6 +104,8 @@ class GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients : public Likel private: const GslBlockMatrix & m_covariance; + + void checkDimConsistency(const V & observations, const GslBlockMatrix & covariance) const; }; } // End namespace QUESO diff --git a/src/stats/inc/GaussianLikelihoodDiagonalCovariance.h b/src/stats/inc/GaussianLikelihoodDiagonalCovariance.h index 438bf2601..da7a6dcbb 100644 --- a/src/stats/inc/GaussianLikelihoodDiagonalCovariance.h +++ b/src/stats/inc/GaussianLikelihoodDiagonalCovariance.h @@ -54,6 +54,32 @@ class GaussianLikelihoodDiagonalCovariance : public LikelihoodBase { const VectorSet & domainSet, const V & observations, const V & covariance); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + GaussianLikelihoodDiagonalCovariance(const char * prefix, + const VectorSet & domainSet, + const V & observations, + const V & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor virtual ~GaussianLikelihoodDiagonalCovariance(); //@} diff --git a/src/stats/inc/GaussianLikelihoodFullCovariance.h b/src/stats/inc/GaussianLikelihoodFullCovariance.h index 1986d505a..00051f312 100644 --- a/src/stats/inc/GaussianLikelihoodFullCovariance.h +++ b/src/stats/inc/GaussianLikelihoodFullCovariance.h @@ -58,6 +58,31 @@ class GaussianLikelihoodFullCovariance : public LikelihoodBase { const VectorSet & domainSet, const V & observations, const M & covariance, double covarianceCoefficient=1.0); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + GaussianLikelihoodFullCovariance(const char * prefix, + const VectorSet & domainSet, const V & observations, + const M & covariance, double covarianceCoefficient, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor virtual ~GaussianLikelihoodFullCovariance(); //@} diff --git a/src/stats/inc/GaussianLikelihoodFullCovarianceRandomCoefficient.h b/src/stats/inc/GaussianLikelihoodFullCovarianceRandomCoefficient.h index 9a71b96f6..0daef6882 100644 --- a/src/stats/inc/GaussianLikelihoodFullCovarianceRandomCoefficient.h +++ b/src/stats/inc/GaussianLikelihoodFullCovarianceRandomCoefficient.h @@ -63,6 +63,31 @@ class GaussianLikelihoodFullCovarianceRandomCoefficient : public LikelihoodBase< const VectorSet & domainSet, const V & observations, const M & covariance); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + GaussianLikelihoodFullCovarianceRandomCoefficient(const char * prefix, + const VectorSet & domainSet, const V & observations, + const M & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor virtual ~GaussianLikelihoodFullCovarianceRandomCoefficient(); //@} diff --git a/src/stats/inc/GaussianLikelihoodScalarCovariance.h b/src/stats/inc/GaussianLikelihoodScalarCovariance.h index ec1c072e1..83d134205 100644 --- a/src/stats/inc/GaussianLikelihoodScalarCovariance.h +++ b/src/stats/inc/GaussianLikelihoodScalarCovariance.h @@ -54,6 +54,31 @@ class GaussianLikelihoodScalarCovariance : public LikelihoodBase { const VectorSet & domainSet, const V & observations, double covariance); + //! Constructor for likelihood that includes marginalization + /*! + * If the likelihood requires marginalization, the user can provide the pdf of the + * marginal parameter(s) and the integration to be used. Additionally, the user will + * be required to have provided an implementation of + * evaluateModel(const V & domainVector, const V & marginalVector, V & modelOutput). + * + * Mathematically, this likelihood evaluation will be + * \f[ \pi(d|m) = \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) \pi(q_i) w_i\f] + * where \f$ N \f$ is the number of quadrature points. However, the PDF for the + * marginal parameter(s) may be such that it is convenient to interpret it as a + * weighting function for Gaussian quadrature. In that case, then, + * \f[ \int \pi(d|m,q) \pi(q)\; dq \approx \sum_{i=1}^{N} \pi(d|m,q_i) w_i \f] + * If this is the case, the user should set the argument marg_pdf_is_weight_func = true. + * If it is set to false, then the former quadrature equation will be used. + * For example, if the marginal parameter(s) pdf is Gaussian, a Gauss-Hermite quadrature + * rule could make sense (GaussianHermite1DQuadrature). + */ + GaussianLikelihoodScalarCovariance(const char * prefix, + const VectorSet & domainSet, const V & observations, + double covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func); + //! Destructor virtual ~GaussianLikelihoodScalarCovariance(); //@} diff --git a/src/stats/src/GaussianLikelihoodBlockDiagonalCovariance.C b/src/stats/src/GaussianLikelihoodBlockDiagonalCovariance.C index 9b8f9f62c..8e032ebba 100644 --- a/src/stats/src/GaussianLikelihoodBlockDiagonalCovariance.C +++ b/src/stats/src/GaussianLikelihoodBlockDiagonalCovariance.C @@ -39,15 +39,21 @@ GaussianLikelihoodBlockDiagonalCovariance::GaussianLikelihoodBlockDiagonal m_covarianceCoefficients(covariance.numBlocks(), 1.0), m_covariance(covariance) { - unsigned int totalDim = 0; - - for (unsigned int i = 0; i < this->m_covariance.numBlocks(); i++) { - totalDim += this->m_covariance.getBlock(i).numRowsLocal(); - } + this->checkDimConsistency(observations,m_covariance); +} - if (totalDim != observations.sizeLocal()) { - queso_error_msg("Covariance matrix not same dimension as observation vector"); - } +template +GaussianLikelihoodBlockDiagonalCovariance::GaussianLikelihoodBlockDiagonalCovariance( + const char * prefix, const VectorSet & domainSet, + const V & observations, const GslBlockMatrix & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func), + m_covarianceCoefficients(covariance.numBlocks(), 1.0), + m_covariance(covariance) +{ + this->checkDimConsistency(observations,m_covariance); } template @@ -105,6 +111,22 @@ GaussianLikelihoodBlockDiagonalCovariance::lnLikelihood(const V & domainVe return -0.5 * norm2_squared; } +template +void +GaussianLikelihoodBlockDiagonalCovariance::checkDimConsistency( + const V & observations, const GslBlockMatrix & covariance) const +{ + unsigned int totalDim = 0; + + for (unsigned int i = 0; i < covariance.numBlocks(); i++) { + totalDim += covariance.getBlock(i).numRowsLocal(); + } + + if (totalDim != observations.sizeLocal()) { + queso_error_msg("Covariance matrix not same dimension as observation vector"); + } +} + } // End namespace QUESO template class QUESO::GaussianLikelihoodBlockDiagonalCovariance; diff --git a/src/stats/src/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.C b/src/stats/src/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.C index a84c6ca0f..8a959ad90 100644 --- a/src/stats/src/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.C +++ b/src/stats/src/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.C @@ -36,15 +36,20 @@ GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::GaussianLikel : LikelihoodBase(prefix, domainSet, observations), m_covariance(covariance) { - unsigned int totalDim = 0; - - for (unsigned int i = 0; i < this->m_covariance.numBlocks(); i++) { - totalDim += this->m_covariance.getBlock(i).numRowsLocal(); - } + this->checkDimConsistency(observations,m_covariance); +} - if (totalDim != observations.sizeLocal()) { - queso_error_msg("Covariance matrix not same dimension as observation vector"); - } +template +GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients( + const char * prefix, const VectorSet & domainSet, + const V & observations, const GslBlockMatrix & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func), + m_covariance(covariance) +{ + this->checkDimConsistency(observations,m_covariance); } template @@ -103,6 +108,22 @@ GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::lnLikelihood( return -0.5 * norm2_squared - cov_norm_factor; } +template +void +GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients::checkDimConsistency( + const V & observations, const GslBlockMatrix & covariance) const +{ + unsigned int totalDim = 0; + + for (unsigned int i = 0; i < covariance.numBlocks(); i++) { + totalDim += covariance.getBlock(i).numRowsLocal(); + } + + if (totalDim != observations.sizeLocal()) { + queso_error_msg("Covariance matrix not same dimension as observation vector"); + } +} + } // End namespace QUESO template class QUESO::GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients; diff --git a/src/stats/src/GaussianLikelihoodDiagonalCovariance.C b/src/stats/src/GaussianLikelihoodDiagonalCovariance.C index 0c94bd75f..c93d31023 100644 --- a/src/stats/src/GaussianLikelihoodDiagonalCovariance.C +++ b/src/stats/src/GaussianLikelihoodDiagonalCovariance.C @@ -43,6 +43,21 @@ GaussianLikelihoodDiagonalCovariance::GaussianLikelihoodDiagonalCovariance } } +template +GaussianLikelihoodDiagonalCovariance::GaussianLikelihoodDiagonalCovariance( + const char * prefix, const VectorSet & domainSet, + const V & observations, const V & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func), + m_covariance(covariance) +{ + if (covariance.sizeLocal() != observations.sizeLocal()) { + queso_error_msg("Covariance matrix not same size as observation vector"); + } +} + template GaussianLikelihoodDiagonalCovariance::~GaussianLikelihoodDiagonalCovariance() { diff --git a/src/stats/src/GaussianLikelihoodFullCovariance.C b/src/stats/src/GaussianLikelihoodFullCovariance.C index 090f4e2b3..76cb70126 100644 --- a/src/stats/src/GaussianLikelihoodFullCovariance.C +++ b/src/stats/src/GaussianLikelihoodFullCovariance.C @@ -44,6 +44,22 @@ GaussianLikelihoodFullCovariance::GaussianLikelihoodFullCovariance( } } +template +GaussianLikelihoodFullCovariance::GaussianLikelihoodFullCovariance(const char * prefix, + const VectorSet & domainSet, const V & observations, + const M & covariance, double covarianceCoefficient, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func), + m_covarianceCoefficient(covarianceCoefficient), + m_covariance(covariance) +{ + if (covariance.numRowsLocal() != observations.sizeLocal()) { + queso_error_msg("Covariance matrix not same size as observation vector"); + } +} + template GaussianLikelihoodFullCovariance::~GaussianLikelihoodFullCovariance() { diff --git a/src/stats/src/GaussianLikelihoodFullCovarianceRandomCoefficient.C b/src/stats/src/GaussianLikelihoodFullCovarianceRandomCoefficient.C index 65481b12e..7b858a535 100644 --- a/src/stats/src/GaussianLikelihoodFullCovarianceRandomCoefficient.C +++ b/src/stats/src/GaussianLikelihoodFullCovarianceRandomCoefficient.C @@ -42,6 +42,21 @@ GaussianLikelihoodFullCovarianceRandomCoefficient::GaussianLikelihoodFullC } } +template +GaussianLikelihoodFullCovarianceRandomCoefficient::GaussianLikelihoodFullCovarianceRandomCoefficient( + const char * prefix, const VectorSet & domainSet, + const V & observations, const M & covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func), + m_covariance(covariance) +{ + if (covariance.numRowsLocal() != observations.sizeLocal()) { + queso_error_msg("Covariance matrix not same size as observation vector"); + } +} + template GaussianLikelihoodFullCovarianceRandomCoefficient::~GaussianLikelihoodFullCovarianceRandomCoefficient() { diff --git a/src/stats/src/GaussianLikelihoodScalarCovariance.C b/src/stats/src/GaussianLikelihoodScalarCovariance.C index c74ab1f93..007af6094 100644 --- a/src/stats/src/GaussianLikelihoodScalarCovariance.C +++ b/src/stats/src/GaussianLikelihoodScalarCovariance.C @@ -40,6 +40,18 @@ GaussianLikelihoodScalarCovariance::GaussianLikelihoodScalarCovariance( { } +template +GaussianLikelihoodScalarCovariance::GaussianLikelihoodScalarCovariance( + const char * prefix, const VectorSet & domainSet, + const V & observations, double covariance, + typename SharedPtr >::Type & marg_param_pdf, + typename SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func), + m_covariance(covariance) +{ +} + template GaussianLikelihoodScalarCovariance::~GaussianLikelihoodScalarCovariance() { From 50c6604a63d158866f879e78c0e8e422aa388acf Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Thu, 9 Feb 2017 09:16:59 -0500 Subject: [PATCH 7/8] Add likelihood marginalization test This is testing the case where we have a uniform RV for the marginal parameters. The likelihood is just f(m) (data =0, effectively uniform pdf) for ease of analytical integration for the test. --- test/Makefile.am | 1 + test/unit/marg_likelihood.C | 227 ++++++++++++++++++++++++++++++++++++ 2 files changed, 228 insertions(+) create mode 100644 test/unit/marg_likelihood.C diff --git a/test/Makefile.am b/test/Makefile.am index 68ff4c4b6..d4e6e6516 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -132,6 +132,7 @@ unit_driver_SOURCES += unit/rng_cxx11.C unit_driver_SOURCES += unit/rng_boost.C unit_driver_SOURCES += unit/basic_pdfs_cxx11.C unit_driver_SOURCES += unit/basic_pdfs_boost.C +unit_driver_SOURCES += unit/marg_likelihood.C test_boxsubset_centroid_SOURCES = test_centroids/test_boxsubset_centroid.C test_concatenation_centroid_SOURCES = test_centroids/test_concatenation_centroid.C diff --git a/test/unit/marg_likelihood.C b/test/unit/marg_likelihood.C new file mode 100644 index 000000000..618eaf977 --- /dev/null +++ b/test/unit/marg_likelihood.C @@ -0,0 +1,227 @@ +//-----------------------------------------------------------------------bl- +//-------------------------------------------------------------------------- +// +// QUESO - a library to support the Quantification of Uncertainty +// for Estimation, Simulation and Optimization +// +// Copyright (C) 2008-2017 The PECOS Development Team +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the Version 2.1 GNU Lesser General +// Public License as published by the Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc. 51 Franklin Street, Fifth Floor, +// Boston, MA 02110-1301 USA +// +//-----------------------------------------------------------------------el- + +#include "config_queso.h" + +#ifdef QUESO_HAVE_CPPUNIT + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +namespace QUESOTesting +{ + template + class TestlingLikelihoodBase : public QUESO::LikelihoodBase + { + public: + + TestlingLikelihoodBase(const char * prefix, + const QUESO::VectorSet & domainSet, + const V & observations, + typename QUESO::SharedPtr >::Type & marg_param_pdf, + typename QUESO::SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : QUESO::LikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func) + {} + + //! Evaluate the exact ln value of the marginalized likelihood + virtual double lnMargLikelihood( const V & domainVector ) const =0; + + void testLikelihoodValue( const V & domainVector, + const QUESO::LikelihoodBase & likelihood, + double tol ) const + { + double computed_value = likelihood.lnValue(domainVector); + + double exact_value = this->lnMargLikelihood(domainVector); + + double rel_error = std::abs( (computed_value-exact_value)/exact_value ); + + CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, rel_error, tol ); + } + + protected: + + virtual double lnLikelihood(const V & /*domainVector*/, + V & modelOutput) const + { + queso_require_equal_to(1,modelOutput.sizeGlobal()); + return std::log(modelOutput[0]); + } + + }; + + template + class LinearModelLikelihood : public TestlingLikelihoodBase + { + public: + + LinearModelLikelihood(const char * prefix, + const QUESO::VectorSet & domainSet, + const V & observations, + typename QUESO::SharedPtr >::Type & marg_param_pdf, + typename QUESO::SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : TestlingLikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func) + {} + + virtual void evaluateModel(const V & domainVector, + const V & marginalVector, + V & modelOutput) const + { + queso_require_equal_to(1,domainVector.sizeGlobal()); + queso_require_equal_to(3,marginalVector.sizeGlobal()); + queso_require_equal_to(1,modelOutput.sizeGlobal()); + + modelOutput[0] = 3.14*domainVector[0] + 1.1*marginalVector[0] + + 2.2*marginalVector[1] + 3.3*marginalVector[2]; + } + + virtual double lnMargLikelihood( const V & domainVector ) const + { + queso_require_equal_to(1,domainVector.sizeGlobal()); + + return std::log(3.14*domainVector[0] + 1.1/2.0 + 2.2/2.0 + 3.3/2.0); + } + + }; + + template + class MarginalLikelihoodTestBase : public CppUnit::TestCase + { + public: + + void setUp() + { + this->init_env(); + } + + void test_linear_func_uniform_marg_space() + { + // Instantiate the parameter space + unsigned int param_dim = 1; + unsigned int marg_dim = 3; + QUESO::VectorSpace param_space( (*this->_env), "param_", param_dim, NULL); + QUESO::VectorSpace marg_space( (*this->_env), "marg_", marg_dim, NULL); + + double param_min_domain_value = 0.0; + double param_max_domain_value = 1.0; + + double marg_min_domain_value = 0.0; + double marg_max_domain_value = 1.0; + + typename QUESO::ScopedPtr::Type param_min_values( param_space.newVector(param_min_domain_value) ); + typename QUESO::ScopedPtr::Type param_max_values( param_space.newVector(param_max_domain_value) ); + + typename QUESO::ScopedPtr::Type marg_min_values( marg_space.newVector(marg_min_domain_value) ); + typename QUESO::ScopedPtr::Type marg_max_values( marg_space.newVector(marg_max_domain_value) ); + + QUESO::BoxSubset param_domain( "param_domain_", param_space, (*param_min_values), (*param_max_values) ); + QUESO::BoxSubset marg_domain( "marg_domain_", marg_space, (*marg_min_values), (*marg_max_values) ); + + typename QUESO::SharedPtr >::Type + marg_param_rv( new QUESO::UniformVectorRV("marg_param_rv_", marg_domain) ); + + const V & data = param_space.zeroVector(); + + unsigned int int_order = 1; + + QUESO::SharedPtr::Type + qrule_1d( new QUESO::UniformLegendre1DQuadrature(marg_min_domain_value, marg_max_domain_value, + int_order, false) ); + + std::vector::Type> all_1d_qrules(marg_dim,qrule_1d); + + typename QUESO::SharedPtr >::Type + marg_integration( new QUESO::TensorProductQuadrature(marg_domain,all_1d_qrules) ); + + LinearModelLikelihood likelihood( "likelihood_test_", + param_domain, + data, + marg_param_rv, + marg_integration, + false ); + + // Test marginalized likelihood over a range of points in the parameter domain + unsigned int n_intervals = 20; + double tol = std::numeric_limits::epsilon()*10; + + this->test_likelihood_values_range( n_intervals, param_min_domain_value, param_max_domain_value, + param_space, tol, likelihood ); + } + + protected: + + QUESO::EnvOptionsValues _options; + + typename QUESO::ScopedPtr::Type _env; + + void init_env() + { + _env.reset( new QUESO::FullEnvironment("","",&_options) ); + } + + void test_likelihood_values_range( unsigned int n_intervals, double param_min_domain_value, + double param_max_domain_value, const QUESO::VectorSpace & param_space, + double tol, const TestlingLikelihoodBase & likelihood ) + { + for( unsigned int i = 0; i < n_intervals; i++ ) + { + double x = param_min_domain_value + i*(param_max_domain_value-param_min_domain_value)/(double)n_intervals; + + typename QUESO::ScopedPtr::Type param_vec( param_space.newVector(x) ); + + likelihood.testLikelihoodValue(*param_vec,likelihood,tol); + } + } + }; + + class MarginalLikelihoodGslTest : public MarginalLikelihoodTestBase + { + public: + CPPUNIT_TEST_SUITE( MarginalLikelihoodGslTest ); + + CPPUNIT_TEST( test_linear_func_uniform_marg_space ); + + CPPUNIT_TEST_SUITE_END(); + }; + + CPPUNIT_TEST_SUITE_REGISTRATION( MarginalLikelihoodGslTest ); + +} // end namespace QUESOTesting + +#endif // QUESO_HAVE_CPPUNIT From e339a52425b0d79d0d7655976b1489b6c9e852d7 Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Thu, 9 Feb 2017 09:54:15 -0500 Subject: [PATCH 8/8] Add another likelihood marginalization test This one is now testing where we treat the marginal parameter pdf as the weighting function in the quadrature rule. Here, we use a zero mean/unit variance Gaussian for the 2D marginal parameter space and then use Gauss-Hermite quadrature. --- test/unit/marg_likelihood.C | 100 ++++++++++++++++++++++++++++++++++++ 1 file changed, 100 insertions(+) diff --git a/test/unit/marg_likelihood.C b/test/unit/marg_likelihood.C index 618eaf977..d67c701f4 100644 --- a/test/unit/marg_likelihood.C +++ b/test/unit/marg_likelihood.C @@ -35,6 +35,7 @@ #include #include #include +#include #include #include @@ -120,6 +121,46 @@ namespace QUESOTesting }; + template + class GaussModelLikelihood : public TestlingLikelihoodBase + { + public: + + GaussModelLikelihood(const char * prefix, + const QUESO::VectorSet & domainSet, + const V & observations, + typename QUESO::SharedPtr >::Type & marg_param_pdf, + typename QUESO::SharedPtr >::Type & marg_integration, + bool marg_pdf_is_weight_func) + : TestlingLikelihoodBase(prefix,domainSet,observations,marg_param_pdf,marg_integration,marg_pdf_is_weight_func) + {} + + virtual void evaluateModel(const V & domainVector, + const V & marginalVector, + V & modelOutput) const + { + queso_require_equal_to(1,domainVector.sizeGlobal()); + queso_require_equal_to(2,marginalVector.sizeGlobal()); + queso_require_equal_to(1,modelOutput.sizeGlobal()); + + double q1 = marginalVector[0]; + double q2 = marginalVector[1]; + double m = domainVector[0]; + + modelOutput[0] = 3.14*m + 1.1*q1 + 2.2*q2*q2; + } + + virtual double lnMargLikelihood( const V & domainVector ) const + { + queso_require_equal_to(1,domainVector.sizeGlobal()); + + double m = domainVector[0]; + + return std::log(3.14*m*M_PI + 2.2*M_PI/2.0); + } + + }; + template class MarginalLikelihoodTestBase : public CppUnit::TestCase { @@ -184,6 +225,64 @@ namespace QUESOTesting param_space, tol, likelihood ); } + void test_linear_func_gaussian_marg_space() + { + // Instantiate the parameter space + unsigned int param_dim = 1; + unsigned int marg_dim = 2; + QUESO::VectorSpace param_space( (*this->_env), "param_", param_dim, NULL); + QUESO::VectorSpace marg_space( (*this->_env), "marg_", marg_dim, NULL); + + double param_min_domain_value = 0.0; + double param_max_domain_value = 1.0; + + double marg_min_domain_value = -INFINITY; + double marg_max_domain_value = INFINITY; + + typename QUESO::ScopedPtr::Type param_min_values( param_space.newVector(param_min_domain_value) ); + typename QUESO::ScopedPtr::Type param_max_values( param_space.newVector(param_max_domain_value) ); + + typename QUESO::ScopedPtr::Type marg_min_values( marg_space.newVector(marg_min_domain_value) ); + typename QUESO::ScopedPtr::Type marg_max_values( marg_space.newVector(marg_max_domain_value) ); + + QUESO::BoxSubset param_domain( "param_domain_", param_space, (*param_min_values), (*param_max_values) ); + QUESO::BoxSubset marg_domain( "marg_domain_", marg_space, (*marg_min_values), (*marg_max_values) ); + + typename QUESO::ScopedPtr::Type var_vec( param_space.newVector(1.0) ); + + // Zero mean, unit variance Gaussian + typename QUESO::SharedPtr >::Type + marg_param_rv( new QUESO::GaussianVectorRV("marg_param_rv_", marg_domain, + marg_space.zeroVector(), + *var_vec) ); + + const V & data = param_space.zeroVector(); + + unsigned int int_order = 1; + + QUESO::SharedPtr::Type + qrule_1d( new QUESO::GaussianHermite1DQuadrature(0.0,1.0,int_order) ); + + std::vector::Type> all_1d_qrules(marg_dim,qrule_1d); + + typename QUESO::SharedPtr >::Type + marg_integration( new QUESO::TensorProductQuadrature(marg_domain,all_1d_qrules) ); + + GaussModelLikelihood likelihood( "likelihood_test_", + param_domain, + data, + marg_param_rv, + marg_integration, + true ); + + // Test marginalized likelihood over a range of points in the parameter domain + unsigned int n_intervals = 20; + double tol = std::numeric_limits::epsilon()*50; + + this->test_likelihood_values_range( n_intervals, param_min_domain_value, param_max_domain_value, + param_space, tol, likelihood ); + } + protected: QUESO::EnvOptionsValues _options; @@ -216,6 +315,7 @@ namespace QUESOTesting CPPUNIT_TEST_SUITE( MarginalLikelihoodGslTest ); CPPUNIT_TEST( test_linear_func_uniform_marg_space ); + CPPUNIT_TEST( test_linear_func_gaussian_marg_space ); CPPUNIT_TEST_SUITE_END(); };