STK_Stat_BivariateRealReal.h

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2011  Serge Iovleff

    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 : Serge.Iovleff@stkpp.org
*/

/*
 * Project:  stkpp::StatistiK::StatDesc
 * Purpose:  Compute elementary statistics for two variables.
 * Author:   Serge Iovleff, serge.iovleff@stkpp.org
 **/

#ifndef STK_STAT_BIVARIATEREALREAL_H
#define STK_STAT_BIVARIATEREALREAL_H

#include "STK_Stat_Bivariate.h"
#include "STK_Stat_UnivariateReal.h"

namespace STK
{
namespace Stat
{

template < class TContainer1D >
class Bivariate<Real, Real, TContainer1D>
{
  public:
    Bivariate( ITContainer1D<Real, TContainer1D> const& X
             , ITContainer1D<Real, TContainer1D> const& Y
             )
             : X_(X.asLeaf())
             , Y_(Y.asLeaf())
             , W_()
             , xStat_(X)
             , yStat_(Y)
             , cov_(0.)
             , ucov_(0.)
             , cor_(0.)
             , ucor_(0.)
             , sumWeights_(0.0)
             , sum2Weights_(0.0)
    {
      compCovariance();
    }

    Bivariate( ITContainer1D<Real, TContainer1D> const& X
             , ITContainer1D<Real, TContainer1D> const& Y
             , TContainer1D const& W
             )
             : X_(X.asLeaf())
             , Y_(Y.asLeaf())
             , W_(W)
             , xStat_(X, W)
             , yStat_(Y, W)
             , cov_(0.)
             , ucov_(0.)
             , cor_(0.)
             , ucor_(0.)
             , sumWeights_(0.0)
             , sum2Weights_(0.0)
    {
      compWeightedCovariance();
    }
    Bivariate( Bivariate<Real, Real, TContainer1D> const& stat)
              : xStat_(stat.xStat_)
              , yStat_(stat.yStat_)
              , cov_(stat.cov_)
              , ucov_(stat.ucov_)
              , cor_(stat.cor_)
              , ucor_(stat.ucor_)
              , sumWeights_(stat.nweight_)
              , sum2Weights_(stat.nweight2_)
    { }
    
    virtual ~Bivariate() { ;}
  
    Bivariate& operator=( Bivariate<Real, Real, TContainer1D> const& stat)
    {
      xStat_ = stat.xStat_;
      yStat_ = stat.yStat_;
      W_     = stat.W_;
      cov_   = stat.cov_;
      ucov_  = stat.ucov_;
      cor_   = stat.cor_;
      ucor_  = stat.ucor_;
      return *this;
    }

    void setData( ITContainer1D<Real, TContainer1D> const& X
                , ITContainer1D<Real, TContainer1D> const& Y
                )
    {
      X_ = X.asLeaf();
      Y_ = Y.asLeaf();
      W_.clear();
      xStat_.setData(X);
      yStat_.setData(Y);
      compCovariance();
    }

    void setData( ITContainer1D<Real, TContainer1D> const& X
                , ITContainer1D<Real, TContainer1D> const& Y
                , TContainer1D const& W
                )
    {
      X_ = X.asLeaf();
      Y_ = Y.asLeaf();
      W_ = W;
      xStat_.setData(X, W);
      yStat_.setData(Y, W);
      compWeightedCovariance();
    }

    inline Univariate<Real, TContainer1D> const& xStat() const {return xStat_;}
    inline Univariate<Real, TContainer1D> const& yStat() const {return yStat_;}
    inline Real const& covariance() const { return cov_;}
    inline Real const& correlation() const { return cor_;}
    inline Real const& unbiasedCovariance() const { return ucov_;}
    inline Real const& unbiasedCorrelation() const { return ucor_;}

  protected:
    TContainer1D X_;
    TContainer1D Y_;
    TContainer1D W_;
    Univariate<Real, TContainer1D>    xStat_;
    Univariate<Real, TContainer1D>    yStat_;
    Real cov_;
    Real ucov_;
    Real cor_;
    Real ucor_;
    Real sumWeights_;
    Real sum2Weights_;

    void compCovariance()
    {
      // get dimensions
      Range com = STK::Range::inf(X_.range(),Y_.range());
      const Integer first = com.first(), last = com.last();
      Integer nobs = com.size();
      // get means
      Real xMean = xStat_.mean(), yMean= yStat_.mean();
      // compute covariance
      Real xsum  = 0.0, ysum = 0.0, xdev, ydev;
      cov_ = 0.;
      for (Integer i=first; i<=last; i++)
      {
        if (Arithmetic<Real>::isFinite(X_[i]) && Arithmetic<Real>::isFinite(Y_[i]))
        {
          xsum += (xdev = X_[i] - xMean); // deviation from the mean
          ysum += (ydev = Y_[i] - yMean); // deviation from the mean
          cov_ += (xdev*ydev);         // squared value
        }
        else nobs--;
      }
      if (nobs > 1)
      {
        ucov_ = (cov_ - (xsum*ysum)/(Real)nobs)/(Real)(nobs -1);
        ucor_ = ucov_/(xStat_.unbiasedStd()*yStat_.unbiasedStd());
        cov_  = (cov_ - (xsum*ysum)/(Real)nobs)/(Real)(nobs);
        cor_  = cov_/(xStat_.std()*yStat_.std());
      }
      else
      {
        if (nobs == 1)
        {
          cov_ = 0.;
          cor_ = 0.;
          ucov_ = Arithmetic<Real>::NA();
          ucor_ = Arithmetic<Real>::NA();
        }
        // no observations
        cov_ = Arithmetic<Real>::NA();
        cor_ = Arithmetic<Real>::NA();
        ucov_ = Arithmetic<Real>::NA();
        ucor_ = Arithmetic<Real>::NA();
      }
    }

    void compWeightedCovariance()
    {
      // get dimensions
      Range com = STK::Range::inf(X_.range(),Y_.range());
      Integer first = com.first(), last = com.last();
      // get means
      Real xMean = xStat_.mean(), yMean= yStat_.mean();
      // compute covariance
      Real xsum  = 0.0, ysum = 0.0, xdev, ydev;
      sumWeights_ = 0.0;
      sum2Weights_ = 0.0;
      // sum
      cov_ = 0.0;
      for (Integer i=first; i<=last; i++)
      {
        if (Arithmetic<Real>::isFinite(X_[i]) && Arithmetic<Real>::isFinite(Y_[i]))
        {
          xsum += (xdev = X_[i] - xMean); // deviation from the mean
          ysum += (ydev = Y_[i] - yMean); // deviation from the mean
          Real Wi = abs((Real)W_[i]);
          cov_ += Wi * (xdev*ydev);         // cross product
          sumWeights_  += Wi;           // sum absolute weights
          sum2Weights_ += Wi * Wi;     // sum squared weights
       }
     }
     // not pathological weights
     if (sumWeights_*sumWeights_ > sum2Weights_)
     {
       ucov_ = (cov_ - xsum*ysum/sumWeights_)/(sumWeights_ - sum2Weights_/sumWeights_);
       ucor_ = ucov_/(xStat_.unbiasedStd()*yStat_.unbiasedStd());
       cov_ = (cov_ - xsum*ysum)/(sumWeights_);
       cor_ = cov_/(xStat_.std()*yStat_.std());
     }
     else
     {
       if (sumWeights_) // if sum of the weights is not 0
       {
         cov_ = 0.;
         cor_ = 0.;
         ucov_ = Arithmetic<Real>::NA();
         ucor_ = Arithmetic<Real>::NA();
       }
       // no weigths
       cov_ = Arithmetic<Real>::NA();
       cor_ = Arithmetic<Real>::NA();
       ucov_ = Arithmetic<Real>::NA();
       ucor_ = Arithmetic<Real>::NA();
     }
   }
};
  
template<class TContainer1D>
Real covarianceWithFixedMean( ITContainer1D<Real, TContainer1D> const& X
                            , ITContainer1D<Real, TContainer1D> const& Y
                            , Real const& xMean
                            , Real const& yMean
                            , bool unbiased = false
                            )
{
    // get dimensions
    Range com = STK::Range::inf(X.range(),Y.range());
    // no samples
    if (com.empty()) { return Arithmetic<Real>::NA();}

    // get dimensions
    const Integer first = com.first(), last = com.last();
    Integer nobs = com.size();
    // compute covariance
    Real xsum  = 0.0, ysum = 0.0, cov  = 0.0, xdev, ydev;
    for (Integer i=first; i<=last; i++)
    {
      if (Arithmetic<Real>::isFinite(X[i]) && Arithmetic<Real>::isFinite(Y[i]))
      {
        xsum += (xdev = X[i] - xMean); // deviation from the mean
        ysum += (ydev = Y[i] - yMean); // deviation from the mean
        cov += (xdev*ydev);         // squared value
      }
      else nobs--;
    }
    // compute the variance
    if (unbiased)
    {
      return (nobs > 1) ? ((cov - (xsum*ysum)/(Real)nobs)/(Real)(nobs -1))
                        : Arithmetic<Real>::NA();

    }
    return (nobs > 0) ? (cov - (xsum*ysum)/(Real)nobs)/(Real)(nobs)
                      : Arithmetic<Real>::NA();
}

template<class TContainer1D>
Real covarianceWithFixedMean( ITContainer1D<Real, TContainer1D> const& X
                            , ITContainer1D<Real, TContainer1D> const& Y
                            , ITContainer1D<Real, TContainer1D> const& W
                            , Real const& xMean
                            , Real const& yMean
                            , bool unbiased = false
                          )
{
    // get dimensions
    Range com = STK::Range::inf(X.range(),Y.range());
    // no samples
    if (com.empty()) { return Arithmetic<Real>::NA();}

    // if the weight are not of the same size, ignore them
    if (!com.isIn(W.range()))
      return covarianceWithFixedMean(X, Y, xMean, yMean, unbiased);

    // get dimensions
    Integer first = com.first(), last = com.last();
    // compute covariance
    Real xsum  = 0.0, ysum = 0.0, xdev, ydev, sumWeights = 0.0, sum2Weights = 0.0, cov = 0.0;
    for (Integer i=first; i<=last; i++)
    {
      if (Arithmetic<Real>::isFinite(X[i]) && Arithmetic<Real>::isFinite(Y[i]))
      {
        xsum += (xdev = X[i] - xMean); // deviation from the mean
        ysum += (ydev = Y[i] - yMean); // deviation from the mean
        Real Wi = abs((Real)W[i]);
        cov += Wi * (xdev*ydev);         // cross product
        sumWeights  += Wi;           // sum absolute weights
        sum2Weights += Wi * Wi;     // sum squared weights
      }
    }
    // compute the variance
    if (unbiased)
    {
      if (sumWeights*sumWeights > sum2Weights)
      {
        return (cov - xsum*ysum/sumWeights)/(sumWeights - sum2Weights/sumWeights);
      }
      else
      {
        if (sumWeights) return 0.;
      }
      return Arithmetic<Real>::NA();
    }
    return (sumWeights) ? (cov - xsum*ysum)/(sumWeights) : Arithmetic<Real>::NA();
}

template<class TContainer1D>
Real covariance( ITContainer1D<Real, TContainer1D> const& X
                , ITContainer1D<Real, TContainer1D> const& Y
                , bool unbiased = false
                )
{
    // get dimensions
    Range com = STK::Range::inf(X.range(),Y.range());
    // no samples
    if (com.empty()) { return Arithmetic<Real>::NA();}
    
    // Compute the mean
    Real xMean = mean<TContainer1D>(X), yMean = mean<TContainer1D>(Y);
    
    // get dimensions
    const Integer first = com.first(), last = com.last();
    Integer nobs = com.size();
    // compute covariance
    Real xsum  = 0.0, ysum = 0.0, cov  = 0.0, xdev, ydev;
    for (Integer i=first; i<=last; i++)
    {
        if (Arithmetic<Real>::isFinite(X[i]) && Arithmetic<Real>::isFinite(Y[i]))
        {
            xsum += (xdev = X[i] - xMean); // deviation from the mean
            ysum += (ydev = Y[i] - yMean); // deviation from the mean
            cov += (xdev*ydev);         // squared value
        }
        else nobs--;
    }
    // compute the variance
    if (unbiased)
    {
        return (nobs > 1) ? ((cov - (xsum*ysum)/(Real)nobs)/(Real)(nobs -1))
        : Arithmetic<Real>::NA();
        
    }
    return (nobs > 0) ? ((cov - (xsum*ysum)/(Real)nobs)/(Real)(nobs))
    : Arithmetic<Real>::NA();
}

template<class TContainer1D>
Real covariance( ITContainer1D<Real, TContainer1D> const& X
                , ITContainer1D<Real, TContainer1D> const& Y
                , ITContainer1D<Real, TContainer1D> const& W
                , bool unbiased = false
                )
{
    // get dimensions
    Range com = STK::Range::inf(X.range(),Y.range());
    // no samples
    if (com.empty()) { return Arithmetic<Real>::NA();}
    
    // Compute the mean
    Real xMean = mean(X, W), yMean= mean(Y, W);
    
    // if the weight are not of the same size, ignore them
    if (!com.isIn(W.range()))
        return covarianceWithFixedMean(X, Y, xMean, yMean, unbiased);
    
    // get dimensions
    Integer first = com.first(), last = com.last();
    // compute covariance
    Real xsum  = 0.0, ysum = 0.0, xdev, ydev, sumWeights = 0.0, sum2Weights = 0.0, cov = 0.0;
    for (Integer i=first; i<=last; i++)
    {
        if (Arithmetic<Real>::isFinite(X[i]) && Arithmetic<Real>::isFinite(Y[i]))
        {
            xsum += (xdev = X[i] - xMean); // deviation from the mean
            ysum += (ydev = Y[i] - yMean); // deviation from the mean
            Real Wi = abs((Real)W[i]);
            cov += Wi * (xdev*ydev);         // cross product
            sumWeights  += Wi;           // sum absolute weights
            sum2Weights += Wi * Wi;     // sum squared weights
        }
    }
    // compute the variance
    if (unbiased)
    {
        if (sumWeights*sumWeights > sum2Weights)
        {
            return (cov - xsum*ysum/sumWeights)/(sumWeights - sum2Weights/sumWeights);
        }
        else
        {
            if (sumWeights) return 0.;
        }
        return Arithmetic<Real>::NA();
    }
    return (sumWeights) ? (cov - xsum*ysum)/(sumWeights) : Arithmetic<Real>::NA();
}


}  // namespace Stat

}  // namespace STK

#endif /*STK_STAT_BIVARIATEREALREAL_H*/