STK_JointGammaModel.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_JOINTGAMMAMODEL_H
#define STK_JOINTGAMMAMODEL_H
 
#include <cmath>
 
#include "STK_IMultiStatModel.h"
#include "STK_IMultiParameters.h"
#include <Arrays/include/STK_Array2DPoint.h>
 
#include <STatistiK/include/STK_Law_Gamma.h>
#include <Analysis/include/STK_Algo_FindZero.h>
#include <Analysis/include/STK_Funct_raw.h>
 
namespace STK
{
struct JointGammaParameters: public IMultiParameters<JointGammaParameters>
{
  JointGammaParameters() : shape_(), scale_(), mean_(), meanLog_(), variance_() {}
  JointGammaParameters( Range const& range)
                      : shape_(range, 1.), scale_(range, 1.)
                      , mean_(range, 1.), meanLog_(range, 0.), variance_(range, 1.)
  {}
  JointGammaParameters( JointGammaParameters const& param)
                      : shape_(param.shape_), scale_(param.scale_)
                      , mean_(param.mean_), meanLog_(param.meanLog_), variance_(param.variance_)
  {}
  ~JointGammaParameters() {}
  inline Real shape(int j) const {return shape_[j];}
  inline Real scale(int j) const {return scale_[j];}
  inline void resizeImpl(Range const& range)
  {
    shape_.resize(range); shape_ = 1.;
    scale_.resize(range); scale_ = 1.;
    mean_.resize(range); mean_ = 1.;
    meanLog_.resize(range); meanLog_ = 0.;
    variance_.resize(range); variance_ = 1.;
  }
  Array2DPoint<Real> shape_;
  Array2DPoint<Real> scale_;
  Array2DPoint<Real> mean_;
  Array2DPoint<Real> meanLog_;
  Array2DPoint<Real> variance_;
};
 
 
template <class Array, class WColVector>
class JointGammaModel: public IMultiStatModel<Array, WColVector, JointGammaParameters >
{
 
  public:
    typedef typename Array::Type Type;
    typedef typename hidden::Traits<Array>::Row RowVector;
    typedef typename hidden::Traits<Array>::Col ColVector;
    typedef IMultiStatModel<Array, WColVector, JointGammaParameters > Base;
    using Base::p_data;
    using Base::p_param;
    JointGammaModel(): Base() {}
    JointGammaModel(Array const& data): Base(data){}
    JointGammaModel(Array const* p_data): Base(p_data){}
    JointGammaModel(JointGammaModel const& model): Base(model) {}
    virtual ~JointGammaModel(){}
    JointGammaModel* clone() const { return new JointGammaModel(*this);}
 
    inline Array2DPoint<Real> const& shape() const {return p_param()->shape_;}
    inline Array2DPoint<Real> const& scale() const {return p_param()->scale_;}
    inline Array2DPoint<Real> const& mean() const {return p_param()->mean_;}
    inline Array2DPoint<Real> const& meanLog() const {return p_param()->meanLog_;}
    inline Array2DPoint<Real> const& variance() const {return p_param()->variance_;}
 
    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::Gamma::lpdf(rowData[j], shape()[j], scale()[j]);}
      return sum;
    }
  protected:
    class dloglikelihood: public IFunction<dloglikelihood >
    {
      public:
        dloglikelihood( Real const& mean, Real const& meanLog)
                      : delta_(meanLog - Real(std::log(mean))) {}
        inline Real fImpl(Real a) const
        { return (delta_ + std::log(a) - Funct::psi_raw(a));}
        inline Real xminImpl() const { return 0;}
      private:
        Real delta_;
    };
    virtual void computeParameters()
    {
      for (int j=p_data()->beginCols(); j < p_data()->endCols(); ++j)
      {
        mean()[j] =  p_data()->col(j).safe().mean();
        meanLog()[j] = p_data()->col(j).safe(1.).log().mean();
        variance()[j] = p_data()->col(j).safe().variance();
        Real x0 = (mean()[j]*mean()[j]) / variance()[j];
        Real x1 = 0.9*x0 +  0.05/(mean()[j] - meanLog()[j]);
        dloglikelihood funct(mean()[j], meanLog()[j]);
        Real a =  Algo::findZero(funct, x0, x1);
        // replace with moment estimate if needed
        if (!Arithmetic<Real>::isFinite(a)) { a =  mean()[j]*mean()[j]/variance()[j];}
        shape()[j] = a;
        scale()[j] = mean()[j]/a;
      }
    }
    virtual void computeParameters( WColVector const& weights)
    {
      for (int j=p_data()->beginCols(); j < p_data()->endCols(); ++j)
      {
        mean()[j] =  p_data()->col(j).safe().wmean(weights);
        meanLog()[j] = p_data()->col(j).safe(1).log().wmean(weights);
        variance()[j] = p_data()->col(j).safe().wvariance(weights);
        Real x0 = (mean()[j]*mean()[j]) / variance()[j];
        Real x1 = 0.9*x0 +  0.05/(mean()[j] - meanLog()[j]);
        dloglikelihood funct(mean()[j], meanLog()[j]);
        Real a =  Algo::findZero(funct, x0, x1);
        // replace with moment estimate if needed
        if (!Arithmetic<Real>::isFinite(a)) { a =  mean()[j]*mean()[j]/variance()[j];}
        shape()[j] = a;
        scale()[j] = mean()[j]/a;
      }
    }
 
  private:
    inline Array2DPoint<Real>& shape() {return p_param()->shape_;}
    inline Array2DPoint<Real>& scale() {return p_param()->scale_;}
    inline Array2DPoint<Real>& mean() {return p_param()->mean_;}
    inline Array2DPoint<Real>& meanLog() {return p_param()->meanLog_;}
    inline Array2DPoint<Real>& variance() {return p_param()->variance_;}
};
 
} // namespace STK
 
#endif /* STK_JOINTGAMMAMODEL_H */