Add Capability to Marginalize Likelihood

src/stats/inc/GaussianLikelihoodBlockDiagonalCovariance.h

@@ -60,6 +60,32 @@ class GaussianLikelihoodBlockDiagonalCovariance : public LikelihoodBase<V, M> {
       const VectorSet<V, M> & 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<V, M> & domainSet,
+                 const V & observations,
+                 const GslBlockMatrix & covariance,
+                 typename SharedPtr<BaseVectorRV<V,M> >::Type & marg_param_pdf,
+                 typename SharedPtr<MultiDQuadratureBase<V,M> >::Type & marg_integration,
+                 bool marg_pdf_is_weight_func);
+
   //! Destructor
   virtual ~GaussianLikelihoodBlockDiagonalCovariance();
   //@}
@@ -70,12 +96,16 @@ class GaussianLikelihoodBlockDiagonalCovariance : public LikelihoodBase<V, M> {
   //! 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<double> m_covarianceCoefficients;
   const GslBlockMatrix & m_covariance;
+
+  void checkDimConsistency(const V & observations, const GslBlockMatrix & covariance) const;
 };
 
 }  // End namespace QUESO

src/stats/inc/GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients.h

@@ -62,24 +62,50 @@ class GaussianLikelihoodBlockDiagonalCovarianceRandomCoefficients : public Likel
       const VectorSet<V, M> & 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<V, M> & domainSet,
+                 const V & observations,
+                 const GslBlockMatrix & covariance,
+                 typename SharedPtr<BaseVectorRV<V,M> >::Type & marg_param_pdf,
+                 typename SharedPtr<MultiDQuadratureBase<V,M> >::Type & marg_integration,
+                 bool marg_pdf_is_weight_func);
+
   //! Destructor
   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;
+
+  void checkDimConsistency(const V & observations, const GslBlockMatrix & covariance) const;
 };
 
 }  // End namespace QUESO

src/stats/inc/GaussianLikelihoodDiagonalCovariance.h

@@ -54,12 +54,40 @@ class GaussianLikelihoodDiagonalCovariance : public LikelihoodBase<V, M> {
       const VectorSet<V, M> & 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<V, M> & domainSet,
+                 const V & observations,
+                 const V & covariance,
+                 typename SharedPtr<BaseVectorRV<V,M> >::Type & marg_param_pdf,
+                 typename SharedPtr<MultiDQuadratureBase<V,M> >::Type & marg_integration,
+                 bool marg_pdf_is_weight_func);
+
   //! Destructor
   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;

src/stats/inc/GaussianLikelihoodFullCovariance.h

@@ -58,12 +58,39 @@ class GaussianLikelihoodFullCovariance : public LikelihoodBase<V, M> {
       const VectorSet<V, M> & 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<V, M> & domainSet, const V & observations,
+      const M & covariance, double covarianceCoefficient,
+      typename SharedPtr<BaseVectorRV<V,M> >::Type & marg_param_pdf,
+      typename SharedPtr<MultiDQuadratureBase<V,M> >::Type & marg_integration,
+      bool marg_pdf_is_weight_func);
+
   //! Destructor
   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;

src/stats/inc/GaussianLikelihoodFullCovarianceRandomCoefficient.h

@@ -63,12 +63,39 @@ class GaussianLikelihoodFullCovarianceRandomCoefficient : public LikelihoodBase<
       const VectorSet<V, M> & 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<V, M> & domainSet, const V & observations,
+      const M & covariance,
+      typename SharedPtr<BaseVectorRV<V,M> >::Type & marg_param_pdf,
+      typename SharedPtr<MultiDQuadratureBase<V,M> >::Type & marg_integration,
+      bool marg_pdf_is_weight_func);
+
   //! Destructor
   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;

src/stats/inc/GaussianLikelihoodScalarCovariance.h

@@ -54,12 +54,39 @@ class GaussianLikelihoodScalarCovariance : public LikelihoodBase<V, M> {
       const VectorSet<V, M> & 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<V, M> & domainSet, const V & observations,
+      double covariance,
+      typename SharedPtr<BaseVectorRV<V,M> >::Type & marg_param_pdf,
+      typename SharedPtr<MultiDQuadratureBase<V,M> >::Type & marg_integration,
+      bool marg_pdf_is_weight_func);
+
   //! Destructor
   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;

src/stats/inc/LikelihoodBase.h

@@ -28,9 +28,16 @@
 #include <vector>
 #include <cmath>
 #include <queso/ScalarFunction.h>
+#include <queso/SharedPtr.h>
 
 namespace QUESO {
 
+template<class V, class M>
+class BaseVectorRV;
+
+template<class V, class M>
+class MultiDQuadratureBase;
+
 class GslVector;
 class GslMatrix;
 
@@ -59,8 +66,34 @@ class LikelihoodBase : public BaseScalarFunction<V, M> {
       const VectorSet<V, M> & 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<V, M> & domainSet,
+                 const V & observations,
+                 typename SharedPtr<BaseVectorRV<V,M> >::Type & marg_param_pdf,
+                 typename SharedPtr<MultiDQuadratureBase<V,M> >::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
@@ -79,7 +112,8 @@ class LikelihoodBase : public BaseScalarFunction<V, M> {
    * 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
@@ -94,13 +128,51 @@ class LikelihoodBase : public BaseScalarFunction<V, M> {
   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
   { 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<const BaseVectorRV<V,M> >::Type m_marg_param_pdf;
+
+  typename SharedPtr<MultiDQuadratureBase<V,M> >::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
+   * 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