STK_IGaussianModel.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::Model
 * 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_IGAUSSIANMODEL_H
#define STK_IGAUSSIANMODEL_H
 
#include "STK_Model_Util.h"
#include "STK_IStatModel.h"
 
#include <Arrays/include/STK_Array2D.h>
#include <Arrays/include/STK_Array2DSquare.h>
#include <Arrays/include/STK_Array2DVector.h>
#include <Arrays/include/STK_Array2DPoint.h>
 
#include <STatistiK/include/STK_Stat_Functors.h>
 
namespace STK
{
template<class> class IGaussianModel;
 
namespace hidden
{
 
template<class Data_>
struct ModelTraits< IGaussianModel<Data_> >
{
  class Void {};
  typedef Data_ Array;
  typedef Data_ Data;
  typedef Void ParamHandler;
};
 
} // namespace hidden
 
template<class Vector>
Real gaussianLnLikelihood(ExprBase<Vector> const& data, Real const& mu, Real const& sigma)
{
  if (sigma)
  {
    Real scale = 0., std = std::sqrt((double)sigma);
    // compute scale
    scale = ((data-mu)/std).maxElt();
    Real sum = 0;
    if (scale)
    { sum = (((data-mu)/std)/scale).norm2();}
    return - ( 0.5*sum*scale*scale
             + data.size() * (std::log((double)std) + Const::_LNSQRT2PI_)
             );
  }
  // 0 variance
  return -STK::Arithmetic<Real>::infinity();
}
 
template <class Array>
class IGaussianModel: public IStatModel< IGaussianModel<Array> >
{
  typedef IStatModel< IGaussianModel<Array> > Base;
  typedef typename hidden::Traits<Array>::Row RowVector;
  typedef typename hidden::Traits<Array>::Col ColVector;
 
  protected:
    IGaussianModel( Array const* p_data): Base(p_data) {}
    IGaussianModel( Array const& data): Base(data) {}
 
  public:
    virtual ~IGaussianModel() {}
    inline RowVector const& mean() const { return mean_;}
 
  protected:
    RowVector mean_;
    inline void compMean() { mean_ = Stat::mean(*(this->p_dataij_));}
    inline void compWeightedMean(ColVector const& weights)
    { mean_ = Stat::mean(*(this->p_data_), weights);}
    virtual void compCovariance() =0;
    virtual void compWeightedCovariance(ColVector const& weights) =0;
};
 
} // namespace STK
 
#endif /* STK_IGAUSSIANMODEL_H */