STK_GaussianModel.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::StatModel
 * created on: 13 août 2011
 * Purpose:  Create a gaussian statistical model.
 * Author:   iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
 *
 **/
 
#ifndef STK_GAUSSIANMODEL_H
#define STK_GAUSSIANMODEL_H
 
#include "STK_IGaussianModel.h"
#include <STatistiK/include/STK_Stat_MultivariateReal.h>
#include <STatistiK/include/STK_MultiLaw_Normal.h>
 
namespace STK
{
 
template <class Array>
class GaussianModel: public IGaussianModel<Array>
{
  public:
    typedef IGaussianModel<Array> Base;
    using Base::p_dataij_;
 
    using Base::mean_;
 
    typedef typename hidden::Traits<Array>::Col ColVector;
    typedef typename hidden::Traits<Array>::Row RowVector;
    GaussianModel( Array const* p_data);
    GaussianModel( Array const& data);
    ~GaussianModel();
    MultiLaw::Normal<RowVector>* const& p_law() const { return p_law_;}
    bool run();
    bool run(ColVector const& weights);
    inline ArraySquareX const& covariance() const { return cov_;}
 
  protected:
    ArraySquareX cov_;
    void compCovariance();
    void compWeightedCovariance(ColVector const& weights);
    MultiLaw::Normal<RowVector>* p_law_;
};
 
/* constructor */
template <class Array>
GaussianModel<Array>::GaussianModel( Array const* p_dataij)
                                 : Base(p_dataij)
                                  , cov_(p_dataij_->cols())
                                  , p_law_(0)
{
  this->setNbFreeParameter(p_dataij_->sizeCols() + (p_dataij_->sizeCols()* (p_dataij_->sizeCols()-1))/2);
}
 
/* constructor */
template <class Array>
GaussianModel<Array>::GaussianModel( Array const& data)
                                  : Base(data)
                                   , cov_(data.cols())
                                   , p_law_(0)
{ setNbFreeParameter(p_dataij_->sizeCols() + (p_dataij_->sizeCols()* (p_dataij_->sizeCols()-1))/2);}
 
/* destructor */
template <class Array>
GaussianModel<Array>::~GaussianModel()
{ if (p_law_) delete p_law_;}
 
/* implementation of the Gaussian statistical model
 * @return @c true if no error occur and @c false otherwise.
 */
template <class Array>
bool GaussianModel<Array>::run()
{
  // compute the mean
  this->compMean();
  // compute the covariance matrix
  compCovariance();
  // create p_law_ (will be deleted in base class)
  // update gaussian law (will be deleted in base class)
  if (!p_law_) p_law_ = new MultiLaw::Normal<RowVector>(mean_, cov_);
  else p_law_->setParameters(mean_, cov_);
  // compute log likelihood of the gaussian law
  this->setLnLikelihood(p_law_->lnLikelihood(*p_dataij_ ));
  // everything ok
  return true;
}
 
/* implementation of the weighted Gaussian statistical model
 * @param weights the weights of the samples
 * @return @c true if no error occur and @c false otherwise.
 */
template <class Array>
bool GaussianModel<Array>::run(ColVector const& weights)
{
  // compute the mean
  this->compWeightedMean(weights);
  // compute the covariance matrix
  compWeightedCovariance(weights);
  // create p_law_ (will be deleted in base class)
  // update gaussian law (will be deleted in base class)
  if (!p_law_) p_law_ = new MultiLaw::Normal<RowVector>(mean_, cov_);
  else p_law_->setParameters(mean_, cov_);
  // compute log likelihood of the gaussian law
  this->setLnLikelihood(p_law_->lnLikelihood(*p_dataij_ ));
  // everything ok
  return true;
}
 
template <class Array>
void GaussianModel<Array>::compCovariance()
{ cov_ = Stat::covariance(*p_dataij_);}
template <class Array>
void GaussianModel<Array>::compWeightedCovariance(ColVector const& weights)
{ cov_ = Stat::covariance(*p_dataij_, weights);}
 
} // namespace STK
 
#endif /* STK_GAUSSIANMODEL_H */