STK_GaussianAAModel.h

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2016  Serge Iovleff, Université Lille 1, Inria
 
    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as
    published by the Free Software Foundation; either version 2 of the
    License, or (at your option) any later version.
 
    This program 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 program; if not, write to the
    Free Software Foundation, Inc.,
    59 Temple Place,
    Suite 330,
    Boston, MA 02111-1307
    USA
 
    Contact : S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
*/
 
/*
 * Project:  stkpp::AAModels
 * Purpose:  Interface base class for AA models.
 * Author:   Serge Iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
 *
 **/
 
#ifndef STK_GAUSSIANAAMODEL_H
#define STK_GAUSSIANAAMODEL_H
 
#include "STK_IAAModel.h"
#include "Arrays/include/STK_Array2DSquare.h"
 
#include <StatModels/include/STK_Model_Util.h>
#include <StatModels/include/STK_IStatModel.h>
#include <STatistiK/include/STK_MultiLaw_Normal.h>
 
namespace STK
{
 
template<class> class GaussianAAModel;
 
namespace hidden
{
 
template<class Data_>
struct ModelTraits< GaussianAAModel<Data_> >
{
    class Void {};
    typedef Data_ Array;
    typedef Data_ Data;
    typedef Void ParamHandler;
};
 
} // namespace hidden
template<class Array>
class GaussianAAModel: public IAAModel<Array>
                      , public IStatModel< GaussianAAModel<Array> >
{
  public:
    typedef IStatModel< GaussianAAModel<Array> > Base;
    using IAAModel<Array>::p_regressor_;
    using IAAModel<Array>::p_reducer_;
    using IAAModel<Array>::p_reduced_;
    using IAAModel<Array>::p_predicted_;
    using IAAModel<Array>::p_residuals_;
    using IAAModel<Array>::dim;
    GaussianAAModel( Array* p_workData);
    GaussianAAModel( Array& workData);
    inline virtual ~GaussianAAModel() {}
    inline Real const& projectedLnLikelihood() const { return projectedLnLikelihood_;}
    inline Real const& residualLnLikelihood() const { return residualLnLikelihood_;}
    inline ArraySquareX const& projectedCovariance() const { return projectedCovariance_;}
    inline ArraySquareX const& residualCovariance() const { return residualCovariance_;}
    inline Real const& residualVariance() const { return residualVariance_;}
    virtual void setWorkData(Array& workData);
    void computeProjectedCovariance();
    void computeModelParameters();
 
  protected:
    void computeNbFreeParameters();
    void computeResidualCovariance();
    void computeProjectedLnLikelihood();
    void computeResidualLnLikelihood();
 
  private:
    ArraySquareX projectedCovariance_;
    ArraySquareX residualCovariance_;
    Real residualVariance_;
    Real projectedLnLikelihood_;
    Real residualLnLikelihood_;
};
 
/* Constructor.
 *  @param p_workData a pointer on the data set to process
 **/
template<class Array>
GaussianAAModel<Array>::GaussianAAModel( Array* p_workData)
                                : IAAModel<Array>(p_workData)
                                , Base(p_workData)
                                , projectedCovariance_()
                                , residualCovariance_()
                                , residualVariance_(0.)
                                , projectedLnLikelihood_(0.)
                                , residualLnLikelihood_(0.)
{ }
 
// constructor
template<class Array>
GaussianAAModel<Array>::GaussianAAModel( Array& workData)
                                : IAAModel<Array>(workData)
                                , Base(workData)
                                , projectedCovariance_()
                                , residualCovariance_()
                                , residualVariance_(0.)
                                , projectedLnLikelihood_(0.)
                                , residualLnLikelihood_(0.)
{ }
 
/* update the container when the data set is modified. **/
template<class Array>
void GaussianAAModel<Array>::setWorkData(Array& workData)
{
  // update data set and flags for the IAAModel part
  IAAModel<Array>::setWorkData(workData);
  // set dimensions to new size for the IStatModel part
  Base::setData(workData);
}
 
/* compute the ln-likelihood of the model */
template<class Array>
void GaussianAAModel<Array>::computeModelParameters()
{
#ifdef STK_AAMODELS_VERBOSE
  stk_cout << _T("GaussianAAModel::computeModelParameters().\n");
#endif
  // compute the number of free parameters
  computeNbFreeParameters();
  // compute the covariance matrix of the residual and the residual variance
  computeResidualCovariance();
  // compute projected lnLikelihood
  computeProjectedLnLikelihood();
  // compute projected lnLikelihood
  computeResidualLnLikelihood();
  // compute complete nLikelihood
  this->setLnLikelihood(projectedLnLikelihood_ + residualLnLikelihood_);
 
#ifdef STK_AAMODELS_VERBOSE
  stk_cout << _T("GaussianAAModel::computeModelParameters() done.\n");
#endif
}
 
 
/* @brief compute the number of free parameter of the model. **/
template<class Array>
void GaussianAAModel<Array>::computeNbFreeParameters()
{
  this->setNbFreeParameter(p_regressor_->nbFreeParameter() + dim() * (dim()+1)/2 + 1);
}
 
/* compute the ln-likelihood of the model */
template<class Array>
void GaussianAAModel<Array>::computeProjectedLnLikelihood()
{
#ifdef STK_AAMODELS_VERBOSE
  stk_cout << _T("GaussianAAModel::computeProjectedLnLikelihood().\n");
#endif
  // range of the column to use
  Range cols = Range(p_reduced_->beginCols(), dim());
  // create a reference with the first columns of the reduced data
  Array reducedData(*p_reduced_, p_reduced_->rows(), cols);
  // create a reference with the first columns of the reduced data
  ArraySquareX reducedCovariance(projectedCovariance(), cols);
  // compute first part of the ln-likelihood
  PointX mean(cols, 0.);
  MultiLaw::Normal<PointX> normalP(mean, reducedCovariance);
  projectedLnLikelihood_ = normalP.lnLikelihood(reducedData);
 
#ifdef STK_AAMODELS_VERBOSE
  stk_cout << _T("GaussianAAModel::computeProjectedLnLikelihood() done.\n");
  stk_cout << _T("projectedLnLikelihood_ = ") << projectedLnLikelihood_ << _T("\n");
#endif
}
 
/* compute the ln-likelihood of the model */
template<class Array>
void GaussianAAModel<Array>::computeResidualLnLikelihood()
{
#ifdef STK_AAMODELS_VERBOSE
  stk_cout << _T("GaussianAAModel::computeResidualsLnLikelihood().\n");
#endif
  // compute constant part and determinant part of the log-likelihood
  residualLnLikelihood_ = ( Const::_LNSQRT2PI_ + 0.5*std::log(residualVariance_ ))
                          * (dim() - this->nbVariable()) * this->nbSample();
  // compute second part of the log-likelihood
  for (int i=p_residuals_->beginRows(); i<p_residuals_->endRows(); i++)
  { residualLnLikelihood_ -= p_residuals_->row(i).norm2()/(2.*residualVariance_);}
 
#ifdef STK_AAMODELS_VERBOSE
  stk_cout << _T("GaussianAAModel::computeResidualsLnLikelihood() done.\n");
  stk_cout << _T("residualLnLikelihood_ = ") << residualLnLikelihood_ << _T("\n");
#endif
}
 
/* compute the covariance matrix of the projected data set. */
template<class Array>
void GaussianAAModel<Array>::computeProjectedCovariance()
{
#ifdef STK_AAMODELS_DEBUG
  if (!p_reduced_)
  { STKRUNTIME_ERROR_NO_ARG(GaussianAAModel::computeProjectedCovariance,p_reduced_ is not set);}
#endif
#ifdef STK_AAMODELS_VERBOSE
  stk_cout << _T("GaussianAAModel::computeProjectedCovariance().\n");
#endif
  projectedCovariance_ = Stat::covariance(*p_reduced_);
#ifdef STK_AAMODELS_VERBOSE
  stk_cout << _T("GaussianAAModel::computeProjectedCovariance() done.\n");
#endif
}
/* compute the covariance matrix of the residuals. */
template<class Array>
void GaussianAAModel<Array>::computeResidualCovariance()
{
#ifdef STK_AAMODELS_DEBUG
  if (!p_residuals_)
  { STKRUNTIME_ERROR_NO_ARG(GaussianAAModel::computeResidualCovariance,p_residuals_ is not set);}
#endif
#ifdef STK_AAMODELS_VERBOSE
  stk_cout << _T("in GaussianAAModel::computeResidualCovariance().\n");
#endif
  residualCovariance_ = Stat::covariance(*p_residuals_);
  residualVariance_ = (residualCovariance_.trace())/Real(this->nbVariable()-dim());
#ifdef STK_AAMODELS_VERBOSE
  stk_cout << _T("GaussianAAModel::computeResidualCovariance() done.\n");
  stk_cout << _T("residualVariance_ = ") << residualVariance_ << _T("\n");
#endif
}
 
 
} // namespace STK
 
#endif //STK_GAUSSIANAAMODEL_H