STK_JointGaussianModel.h

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2015  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: 22 juil. 2011
 * Purpose: define the class IUnivStatModel.
 * Author:   iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
 *
 **/
 
#ifndef STK_JOINTGAUSSIANMODEL_H
#define STK_JOINTGAUSSIANMODEL_H
 
#include <cmath>
 
#include "STK_IMultiStatModel.h"
#include "STatistiK/include/STK_Law_Normal.h"
#include "STatistiK/include/STK_Stat_Functors.h"
 
namespace STK
{
 
struct JointGaussianParameters: public IMultiParameters<JointGaussianParameters>
{
  public:
    JointGaussianParameters() : mu_(), sigma_() {}
    JointGaussianParameters(Range const& range) : mu_(range, 0.), sigma_(range, 1.) {}
    JointGaussianParameters( JointGaussianParameters const& param)
                            : mu_(param.mu_)
                            , sigma_(param.sigma_)
    {}
    ~JointGaussianParameters() {}
    inline Array2DPoint<Real> const& mu() const { return mu_;}
    inline Array2DPoint<Real> const& sigma() const { return sigma_;}
    inline Real const mu(int const& j) const { return mu_[j];}
    inline Real const sigma(int const& j) const { return sigma_[j];}
    inline void setMu(int const& j, Real const& mu) { mu_[j] = mu;}
    inline void setSigma(int const& j, Real const& sigma) { sigma_[j] = sigma;}
    inline void resizeImpl(Range const& range)
    { mu_.resize(range); mu_ = 0.;
      sigma_.resize(range); sigma_ = 1.;
    }
 
  protected:
    Array2DPoint<Real> mu_;
    Array2DPoint<Real> sigma_;
};
 
template <class Array, class WColVector = CVectorX>
class JointGaussianModel: public IMultiStatModel<Array, WColVector, JointGaussianParameters >
{
  public:
    typedef typename Array::Type Type;
    typedef typename hidden::Traits<Array>::Row RowVector;
    typedef typename hidden::Traits<Array>::Col ColVector;
    typedef IMultiStatModel<Array, WColVector, JointGaussianParameters > Base;
    using Base::p_data;
    using Base::p_param;
    JointGaussianModel(): Base() {}
    JointGaussianModel(Array const& data): Base(data) {}
    JointGaussianModel(Array const* p_data): Base(p_data) {}
    JointGaussianModel(JointGaussianModel const& model): Base(model) {}
    virtual ~JointGaussianModel(){}
    JointGaussianModel* clone() const { return new JointGaussianModel(*this);}
 
    virtual int computeNbFreeParameters() const
    { return 2*p_data()->sizeCols();}
    virtual Real computeLnLikelihood( RowVector const& rowData) const
    {
      Real sum =0.;
      for (Integer j= rowData.begin(); j < rowData.end(); ++j)
      { sum += Law::Normal::lpdf(rowData[j], p_param()->mu(j), p_param()->sigma(j));}
      return sum;
    }
  protected:
    virtual void computeParameters()
    {
      for (int j=p_data()->beginCols(); j < p_data()->endCols(); ++j)
      {
        p_param()->setMu(j, p_data()->col(j).safe().mean());
        p_param()->setSigma(j,std::sqrt(Stat::varianceWithFixedMean(p_data()->col(j), p_param()->mu(j), true)));
      }
    }
    virtual void computeParameters( WColVector const& weights)
    {
      for (int j=p_data()->beginCols(); j < p_data()->endCols(); ++j)
      {
        p_param()->setMu(j, p_data()->col(j).safe().wmean(weights));
        p_param()->setSigma(j,std::sqrt(Stat::varianceWithFixedMean(p_data()->col(j), weights, p_param()->mu(j), true)));
      }
    }
};
 
} // namespace STK
 
#endif /* STK_JOINTGAUSSIANMODEL_H */