STK_ModelDiagGaussian_muj_sj.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: 22 juil. 2011
 * Purpose: define the class IUnivStatModel.
 * Author:   iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
 *
 **/
 
#ifndef STK_GAUSSIAN_MUJ_SJMODEL_H
#define STK_GAUSSIAN_MUJ_SJMODEL_H
 
#include <cmath>
 
#include "STK_Model_Util.h"
#include "STK_IMultiStatModel.h"
#include <STatistiK/include/STK_Law_Normal.h>
#include <STatistiK/include/STK_Stat_Functors.h>
 
namespace STK
{
 
// forward declaration
template <class Data_, class WColVector = CVectorX> class ModelDiagGaussian_muj_sj;
struct DiagGaussian_muj_sjParameters;
 
namespace hidden
{
template<class Data_, class WColVector_>
struct StatModelTraits< ModelDiagGaussian_muj_sj<Data_, WColVector_> >
{
  typedef DataBridge<Data_> Data;
  typedef WColVector_ WColVector;
  typedef typename Traits<Data_>::Type Type;
  typedef DiagGaussian_muj_sjParameters Parameters;
};
 
} // hidden
 
struct DiagGaussian_muj_sjParameters
{
    DiagGaussian_muj_sjParameters(): mu_(), sigma_() {}
    DiagGaussian_muj_sjParameters(Range const& range): mu_(range, 0.), sigma_(range, 1.) {}
    DiagGaussian_muj_sjParameters( DiagGaussian_muj_sjParameters const& param)
                                 : mu_(param.mu_)
                                 , sigma_(param.sigma_)
    {}
    ~DiagGaussian_muj_sjParameters() {}
 
    inline CPointX const& mu() const { return mu_;}
    inline CPointX const& sigma() const { return sigma_;}
 
    inline Real const mu(int j) const { return mu_[j];}
    inline Real const sigma(int j) const { return sigma_[j];}
 
    inline void resize(Range const& range)
    {
      if (range != range_)
      {
        mu_.resize(range); sigma_.resize(range);
        range_ = range;
      }
    }
    CPointX mu_;
    CPointX sigma_;
    Range range_;
};
 
template <class Data_, class WColVector_>
class ModelDiagGaussian_muj_sj: public IMultiStatModel< ModelDiagGaussian_muj_sj<Data_, WColVector_> >
{
  public:
    typedef DataBridge<Data_> Data;
    typedef typename hidden::Traits<Data_>::Row RowVector;
    typedef WColVector_ WColVector;
    typedef typename hidden::Traits<Data_>::Type Type;
    typedef DiagGaussian_muj_sjParameters Parameters;
    typedef IMultiStatModel< ModelDiagGaussian_muj_sj<Data_, WColVector_> > Base;
    using Base::p_data;
    using Base::param;
 
    ModelDiagGaussian_muj_sj(): Base() {}
    ModelDiagGaussian_muj_sj(Data const& data): Base(data) {}
    ModelDiagGaussian_muj_sj(Data const* p_data): Base(p_data) {}
    ModelDiagGaussian_muj_sj(ModelDiagGaussian_muj_sj const& model): Base(model) {}
    ~ModelDiagGaussian_muj_sj() {}
 
    inline CPointX const& mean() const {return param().mu();}
    inline CPointX const& sigma() const {return param().sigma();}
 
    inline int computeNbFreeParameters() const { return 2*p_data()->dataij().sizeCols();}
    Real computeLnLikelihood( RowVector const& rowData) const;
    void computeParameters();
    void computeParameters( WColVector const& weights);
    void writeParametersImpl(ostream& os) const;
};
 
/* compute the log Likelihood of an observation. */
template<class Data_, class WColVector_>
Real ModelDiagGaussian_muj_sj<Data_, WColVector_>::computeLnLikelihood( RowVector const& rowData) const
{
  Real sum =0.;
  for (Integer j= rowData.begin(); j < rowData.end(); ++j)
  { sum += Law::Normal::lpdf(rowData[j], param().mu(j), param().sigma(j));}
  return sum;
}
 
/* compute the parameters */
template<class Data_, class WColVector_>
void ModelDiagGaussian_muj_sj<Data_, WColVector_>::computeParameters()
{
  for (int j=p_data()->dataij().beginCols(); j < p_data()->dataij().endCols(); ++j)
  {
    param().mu_[j] = p_data()->dataij().col(j).meanSafe();
    param().sigma_[j] = std::sqrt(p_data()->dataij().col(j).varianceSafe(param().mu(j)));
  }
}
/* compute the weighted parameters */
template<class Data_, class WColVector_>
void ModelDiagGaussian_muj_sj<Data_, WColVector_>::computeParameters( WColVector const& weights)
{
  for (int j=p_data()->dataij().beginCols(); j < p_data()->dataij().endCols(); ++j)
  {
    param().mu_[j] = p_data()->dataij().col(j).wmeanSafe(weights);
    param().sigma_[j] = std::sqrt(p_data()->dataij().col(j).wvarianceSafe(param().mu(j), weights));
  }
}
 
/* Write the parameters on the output stream os */
template<class Data_, class WColVector_>
void ModelDiagGaussian_muj_sj<Data_, WColVector_>::writeParametersImpl(ostream& os) const
{
  os << _T("mean = ") << mean();
  os << _T("sigma = ")<< sigma();
}
 
} // namespace STK
 
#endif /* STK_GAUSSIAN_MUJ_SJMODEL_H */