STK_Stat_UnivariateReal.h

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2007  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 1D statistics for all variables.
 * Author:   Serge Iovleff, serge.iovleff@stkpp.org
 **/

#ifndef STK_STAT_UNIVARIATEREAL_H
#define STK_STAT_UNIVARIATEREAL_H

#include "../../Sdk/include/STK_ITContainer1D.h"
#include "../../DManager/include/STK_HeapSort.h"

#include "STK_Stat_Univariate.h"

namespace STK
{
namespace Stat
{

template < class TContainer1D >
class Univariate<Real, TContainer1D>
{
  public:
    Univariate( ITContainer1D<Real, TContainer1D> const& V, bool sorted = false)
              : nbSamples_(V.size())
              , nbAvailable_(V.size())
              , nbMiss_(0)
              , V_(V.asLeaf())
              , W_()
              , weighted_(false)
              , sorted_(sorted)
              , comporder_(false)
              , compstat_(false)
              , sumweights_(0.)
              , sum2weights_(0.)
              , mean_(Arithmetic<Real>::NA())
              , median_(Arithmetic<Real>::NA())
              , var_(Arithmetic<Real>::NA())
              , uvar_(Arithmetic<Real>::NA())
              , std_(Arithmetic<Real>::NA())
              , ustd_(Arithmetic<Real>::NA())
              , mad_(Arithmetic<Real>::NA())
              , kurtosis_(Arithmetic<Real>::NA())
              , skewness_(Arithmetic<Real>::NA())
              , quartiles_(3, Arithmetic<Real>::NA())
              , deciles_(9, Arithmetic<Real>::NA())
              , viceciles_(19, Arithmetic<Real>::NA())
              , percentiles_(99, Arithmetic<Real>::NA())
    {
      initializeVariable();
      compOrderStatistics();
      compStatistics();
    }

    Univariate( ITContainer1D<Real, TContainer1D> const& V
              , ITContainer1D<Real, TContainer1D> const& W
              , bool sorted = false
              )
              : nbSamples_(V.size())
              , nbAvailable_(V.size())
              , nbMiss_(0)
              , V_(V.asLeaf())
              , W_(W.asLeaf())
              , weighted_(true)
              , sorted_(sorted)
              , comporder_(false)
              , compstat_(false)
              , sumweights_(0.0)
              , sum2weights_(0.0)
              , mean_(Arithmetic<Real>::NA())
              , median_(Arithmetic<Real>::NA())
              , var_(Arithmetic<Real>::NA())
              , uvar_(Arithmetic<Real>::NA())
              , std_(Arithmetic<Real>::NA())
              , ustd_(Arithmetic<Real>::NA())
              , mad_(Arithmetic<Real>::NA())
              , kurtosis_(Arithmetic<Real>::NA())
              , skewness_(Arithmetic<Real>::NA())
              , quartiles_(3)
              , deciles_(9)
              , viceciles_(19)
              , percentiles_(99)
    {
      // check weights size
      if ((V_.range() != W_.range()))
      {  throw runtime_error("Univariate<Real, TContainer1D>::Univariate(V, W, sorted) "
                                  "V and W have not the same range");
      }
      initializeVariableAndWeights();
      compOrderStatistics();
      compWeightedStatistics();
    }

    Univariate( const Univariate& stat)
              : nbSamples_(stat.nbSamples_)
              , nbAvailable_(stat.nbAvailable_)
              , nbMiss_(stat.nbMiss_)
              , V_(stat.V_)
              , W_(stat.W_)
              , weighted_(stat.weighted_)
              , sorted_(stat.sorted_)
              , comporder_(stat.comporder_)
              , compstat_(stat.compstat_)
              , sumweights_(stat.sumweights_)
              , sum2weights_(stat.sum2weights_)
              , min_(stat.min_)
              , max_(stat.max_)
              , amax_(stat.amax_)
              , mean_(stat.mean_)
              , median_(stat.median_)
              , var_(stat.var_)
              , uvar_(stat.uvar_)
              , std_(stat.std_)
              , ustd_(stat.ustd_)
              , mad_(stat.mad_)
              , kurtosis_(stat.kurtosis_)
              , skewness_(stat.skewness_)
              , quartiles_(stat.quartiles_)
              , deciles_(stat.deciles_)
              , viceciles_(stat.viceciles_)
              , percentiles_(stat.percentiles_)
    { ;}

    ~Univariate() { ;}

    Univariate& operator=( const Univariate& stat)
    {
      nbSamples_ = stat.nbSamples_;
      nbAvailable_ = stat.nbAvailable_;
      nbMiss_ = stat.nbMiss_;
      V_           = stat.V_;
      W_           = stat.W_;
      weighted_    = stat.weighted_;
      sorted_      = stat.sorted_;
      comporder_   = stat.comporder_;
      compstat_    = stat.compstat_;
      sumweights_  = stat.sumweights_;
      sum2weights_ = stat.sum2weights_;
      min_         = stat.min_;
      max_         = stat.max_;
      amax_        = stat.amax_;
      mean_        = stat.mean_;
      median_      = stat.median_;
      var_         = stat.var_;
      uvar_        = stat.uvar_;
      std_         = stat.std_;
      ustd_        = stat.ustd_;
      mad_         = stat.mad_;
      kurtosis_    = stat.kurtosis_;
      skewness_    = stat.skewness_;
      quartiles_   = stat.quartiles_;
      deciles_     = stat.deciles_;
      viceciles_   = stat.viceciles_;
      percentiles_ = stat.percentiles_;

      return *this;
    }

    void setData( ITContainer1D<Real, TContainer1D> const& V, bool sorted = false)
    {
      // initialize dimensions
      nbSamples_    = V.size();
      nbAvailable_  = V.size();
      nbMiss_ = 0;
      // set variable and weights
      V_ = V.asLeaf();
      W_.clear();  // no weights
      // initialize data
      weighted_    = false;
      sorted_      = sorted;
      comporder_   = false;
      compstat_    = false;
      sumweights_     = 0.0;
      sum2weights_    = 0.0;
      mean_        = Arithmetic<Real>::NA();
      median_      = Arithmetic<Real>::NA();
      var_         = Arithmetic<Real>::NA();
      uvar_        = Arithmetic<Real>::NA();
      std_         = Arithmetic<Real>::NA();
      ustd_        = Arithmetic<Real>::NA();
      mad_         = Arithmetic<Real>::NA();
      kurtosis_    = Arithmetic<Real>::NA();
      skewness_    = Arithmetic<Real>::NA();
      quartiles_   = Arithmetic<Real>::NA();
      deciles_     = Arithmetic<Real>::NA();
      viceciles_   = Arithmetic<Real>::NA();
      percentiles_ = Arithmetic<Real>::NA();

      initializeVariable();
      compOrderStatistics();
      compStatistics();
    }

    void setData( ITContainer1D<Real, TContainer1D> const& V
                , ITContainer1D<Real, TContainer1D> const& W
                , bool sorted = false
                )
    {
      // check weights size
      if ((V.range() != W.range()))
      {  throw runtime_error("Univariate<Real>::setData(V, W, sorted) "
                                  "V and W have not the same range");
      }
      // initialize dimensions
      nbSamples_     = V.size();
      nbAvailable_  = V.size();
      nbMiss_ = 0;
      // set variable and weights
      V_ = V.asLeaf();
      W_ = W.asLeaf();
      // initialize data
      weighted_    = true;
      sorted_      = sorted;
      comporder_   = false;
      compstat_    = false;
      sumweights_  = 0.0;
      sum2weights_ = 0.0;
      mean_        = Arithmetic<Real>::NA();
      median_      = Arithmetic<Real>::NA();
      var_         = Arithmetic<Real>::NA();
      uvar_        = Arithmetic<Real>::NA();
      std_         = Arithmetic<Real>::NA();
      ustd_        = Arithmetic<Real>::NA();
      mad_         = Arithmetic<Real>::NA();
      kurtosis_    = Arithmetic<Real>::NA();
      skewness_    = Arithmetic<Real>::NA();
      quartiles_   = Arithmetic<Real>::NA();
      deciles_     = Arithmetic<Real>::NA();
      viceciles_   = Arithmetic<Real>::NA();
      percentiles_ = Arithmetic<Real>::NA();

      initializeVariableAndWeights();
      compOrderStatistics();
      compWeightedStatistics();
    }

    inline const Integer nbSamples() const {return nbSamples_;}
    inline const Integer nbAvailableSamples() const {return nbAvailable_;}
    inline const Integer nbMissingSamples() const {return nbMiss_;}

    inline const Real min() const {return min_;}
    inline const Real max() const {return max_;}
    inline const Real aMax() const {return amax_;}

    inline const Real mean() const {return mean_;}
    inline const Real median() const {return median_;}

    inline const Real variance() const {return var_;}
    inline const Real unbiasedVariance() const {return uvar_;}
    inline const Real std() const {return std_;}
    inline const Real unbiasedStd() const {return ustd_;}
    inline const Real mad() const { return( mad_);}

    inline const Real kurtosis() const {return kurtosis_;}
    inline const Real skewness() const {return skewness_;}

    inline TContainer1D const& quartiles() const {return quartiles_;}
    inline TContainer1D const& deciles() const {return deciles_;}
    inline TContainer1D const& viceciles() const { return viceciles_;}
    inline TContainer1D const& percentiles() const { return percentiles_;}

  private:
    void initializeVariable()
    {
      // initialize the extreme values
      min_ =  Arithmetic<Real>::max();
      max_ = -Arithmetic<Real>::max();
      // discard missing values
      for (Integer i=V_.last(); i>=V_.first(); i--)
      {
        if ( Arithmetic<Real>::isNA(V_[i]) )
        {
          // if NA
          nbAvailable_--;         // decrease nbAvailable_
          nbMiss_++;        // increase nbMiss_
          V_.erase(i);    // Delete the current row
        }
        else
        {
          min_ = STK::min(V_[i], min_);    // update min_
          max_ = STK::max(V_[i], max_);    // update max_
        }
      }
      // all weights are 1/n.
      sumweights_ = 1.0;
      sum2weights_= 1./(Real)nbAvailable_;

      // no samples
      if ((nbAvailable_ == 0))
      {
        min_  = Arithmetic<Real>::NA();
        max_  = Arithmetic<Real>::NA();
        amax_ = Arithmetic<Real>::NA();
      }
      else // one or more samples : compute amax
        amax_ = STK::max(abs(min_), abs(max_));
    }

    void initializeVariableAndWeights()
    {
      // initialize the extreme values
      min_ =  Arithmetic<Real>::max();
      max_ = -Arithmetic<Real>::max();

      // discard missing values or values without weights
      for (Integer i=V_.last(); i>=V_.first(); i--)
      {
        if ( Arithmetic<Real>::isNA(V_[i])||Arithmetic<Real>::isNA(W_[i]))
        {
          V_.erase(i); // Delete the current row
          W_.erase(i); // for the weights too
          nbAvailable_--;         // decrease number of samples
          nbMiss_++;        // increase number of missing values
        }
        else
        {
          min_       = STK::min(V_[i], min_);    // update min_
          max_       = STK::max(V_[i], max_);    // update max_
          sumweights_  += (W_[i] = abs((Real)W_[i]));// sum absolute weights
          sum2weights_ += W_[i] * W_[i];     // sum squared weights
        }
      }
      // no samples
      if ((nbAvailable_ == 0))
      {
        min_  = Arithmetic<Real>::NA();
        max_  = Arithmetic<Real>::NA();
        amax_ = Arithmetic<Real>::NA();
      }
      else // one or more samples : compute amax
        amax_ = STK::max(abs(min_), abs(max_));
    }

    void compStatistics()
    {
      // If there is no samples
      if (nbAvailable_ == 0)
      {
        mean_     = Arithmetic<Real>::NA();
        var_      = Arithmetic<Real>::NA();
        uvar_     = Arithmetic<Real>::NA();
        std_      = Arithmetic<Real>::NA();
        ustd_     = Arithmetic<Real>::NA();
        mad_      = Arithmetic<Real>::NA();
        kurtosis_ = Arithmetic<Real>::NA();
        skewness_ = Arithmetic<Real>::NA();
        return;
      }

      // One observation
      if (nbAvailable_ == 1)
      {
        mean_     = V_.front();
        var_      = Arithmetic<Real>::NA();
        std_      = 0.0;
        ustd_     = Arithmetic<Real>::NA();
        mad_      = 0.0;
        kurtosis_ = 0.0;
        skewness_ = 0.0;
        return;
      }
      // get indexes
      const Integer first = V_.first(), last = V_.last();
      // get nbAvailable observation in Real
      const Real n = (Real)nbAvailable_;
      // first pass : compute the mean
      // scale the samples for preventing overflows
      mean_ = 0.;
      {
        if (amax_) // if the absolute maximal value is greater than 0
        {
          // sum samples with scaling
          for (Integer i=first; i<=last; i++)
            mean_ += V_[i]/amax_;
          // divide by the number of available observation
          mean_ /= n;
          // unscale
          mean_ *= amax_;
        }
      }
      // second pass : compute the variance, skewness and kurtosis
      // initialize values
      Real sum = 0.0, dev1, dev2;
      var_  = 0.0;
      mad_  = 0.0;
      kurtosis_ = 0.0;
      skewness_ = 0.0;
      for (Integer i=1; i<=nbAvailable_; i++)
      {
        sum       += (dev1 = V_[i]-mean_); // deviation from the mean
        var_      += (dev2 = dev1*dev1);   // squared deviation
        mad_      += abs( dev1 );          // absolute deviation
        skewness_ += dev2 * dev1;          // cubic deviation
        kurtosis_ += dev2 * dev2;          // squared squared deviation
      }
      mad_ *= sum2weights_;
      uvar_ = (var_ - sum*sum/n)/(n - 1.);
      var_  = (var_ - sum*sum/n)/n;
      std_  = sqrt((double)var_);
      ustd_ = sqrt((double)uvar_);
     // if there is variance
      if (var_)
      {
        skewness_ /= n;
        kurtosis_ /= n;

        skewness_ /= (var_ * std_);
        kurtosis_ /= (var_ * var_);
        kurtosis_ -= 3.0;
      }
      else
      {
        skewness_ = 0.0;
        kurtosis_ = 0.0;
      }
      // Statistics are computed
      compstat_ = true;
    }

    void compWeightedStatistics()
    {
      // If there is no samples
      if (nbAvailable_ == 0)
      {
        mean_     = Arithmetic<Real>::NA();
        var_      = Arithmetic<Real>::NA();
        uvar_     = Arithmetic<Real>::NA();
        std_      = Arithmetic<Real>::NA();
        ustd_     = Arithmetic<Real>::NA();
        mad_      = Arithmetic<Real>::NA();
        kurtosis_ = Arithmetic<Real>::NA();
        skewness_ = Arithmetic<Real>::NA();
        return;
      }

      // One observation or pathological weights
      if ((nbAvailable_ == 1)||(sumweights_*sumweights_ <= sum2weights_))
      {
        mean_     = V_.front();
        var_      = 0.;
        uvar_     = Arithmetic<Real>::NA();
        std_      = 0.;
        ustd_     = Arithmetic<Real>::NA();
        mad_      = 0.;
        kurtosis_ = 0.;
        skewness_ = 0.;
        return;
      }
      // get indexes
      const Integer first = V_.first(), last = V_.last();
      // first pass : get the mean
      // scale the samples for preventing overflows
      mean_ = 0.;
      if (amax_) // if the maximal value is greater than 0
       {
        // sum samples
        for (Integer i=first; i<=last; i++)
          mean_ += (W_[i]*V_[i])/amax_; // compute the mean with scaling
        mean_ /= sumweights_;              // weight the sum
        mean_ *= amax_;                 // and unscale
      }
      // second pass : compute the variance, skewness and kurtosis
      // initialize some values
      Real sum = 0.0, dev1, dev2;
      var_  = 0.0;
      mad_  = 0.0;
      kurtosis_ = 0.0;
      skewness_ = 0.0;
      for (Integer i=1; i<=nbAvailable_; i++)
      {
        Real weight = W_[i];
        sum       += weight * (dev1 = V_[i]-mean_); // deviation from the mean
        var_      += weight * (dev2 = dev1*dev1);   // squared value
        mad_      += weight * abs(dev1);       // absolute deviation from mean
        skewness_ += weight * dev2 * dev1;
        kurtosis_ += weight * dev2 * dev2;
      }
      mad_ /= sumweights_;
      uvar_ = (var_ - sum*sum/sumweights_)/(sumweights_ - sum2weights_/sumweights_);
      var_  = (var_ - sum*sum)/(sumweights_);
      std_  = sqrt((double)var_);
      ustd_ = sqrt((double)uvar_);
      // if there is variance
      if (var_)
      {
        skewness_ /= (sumweights_ * var_ * std_);
        kurtosis_ /= (sumweights_ * var_ * var_);
        kurtosis_ -= 3.0;
      }
      else
      {
        skewness_ = 0.0;
        kurtosis_ = 0.0;
      }
      // Stat are computed
      compstat_ = true;
    }

    void compOrderStatistics()
    {
      // no samples
      if (nbAvailable_ == 0)
      {
        median_      = Arithmetic<Real>::NA();
        quartiles_   = Arithmetic<Real>::NA();
        deciles_     = Arithmetic<Real>::NA();
        viceciles_   = Arithmetic<Real>::NA();
        percentiles_ = Arithmetic<Real>::NA();
        return;
      }

      // One observation
      if (nbAvailable_ == 1)
      {
        median_      = max_;
        quartiles_   = max_;
        deciles_     = max_;
        viceciles_   = max_;
        percentiles_ = max_;
        return;
      }
      // sort values
      if (!sorted_)
      {
        // if the Variable is not weighted, we can sort it directly
        if (!weighted_) heapSort(V_);
        else // otherwise we have to sort V_and W_ using indirection
        {
          Array1D<Integer> I; // auxiliary Array for indirection
          heapSort(I, V_);
          applySort(V_, I);
          applySort(W_, I);
        }
        sorted_ = true;
      }
      // Find the Quantiles of the distribution
      compQuantiles(quartiles_);
      compQuantiles(deciles_);
      compQuantiles(viceciles_);
      compQuantiles(percentiles_);
      // get the median
      median_ = quartiles_[2];
    }

    void compQuantiles(TContainer1D& T)
    {
#ifdef STK_DEBUG
       if (!sorted_)
      {  throw runtime_error("Univariate<Real>::comptiles(T) "
                                  "V_ is not sorted");
      }
#endif
      // in case of
      T.shift(1);
      // number of quantiles
      Integer  nt = T.last(), shift = V_.first()-1;
      Real n1 = Real(nbAvailable_+1), nt1 = Real(nt+1);

      for (Integer k=1; k<=nt; k++)
      {
        // find index of the k-th quantile
        Real    find  = Real(k)*n1/nt1;  // compute the index in Real
        Integer tind  = STK::max(Integer(1), Integer(find));  // in Integer
        Integer tind1 = STK::min(nbAvailable_, tind+1);              // next
        Real    aux   = find - Real(tind); // lower ponderation

        // nbAvailable_+1 not perfectly divisible ? weighting...
        aux ? T[k] = aux * V_[shift+tind] + (1.0-aux)*V_[shift+tind1]
            : T[k] = V_[shift+tind];
      }
    }

  protected:
    // dimensions
    Integer    nbSamples_;    
    Integer    nbAvailable_; 
    Integer    nbMiss_;      

    // containers
    TContainer1D V_;      
    TContainer1D W_;            

    // Some flag about the internal state of the object
    bool weighted_;       
    bool sorted_;         
    bool comporder_;      
    bool compstat_;       

    // statistics
    Real sumweights_;        
    Real sum2weights_;       
    Real min_;            
    Real max_;            
    Real amax_;           

    Real mean_;           
    Real median_;         
    Real var_;            
    Real uvar_;           
    Real std_;            
    Real ustd_;           
    Real mad_;            
    Real kurtosis_;       
    Real skewness_;       

    TContainer1D quartiles_;    
    TContainer1D deciles_;      
    TContainer1D viceciles_;    
    TContainer1D percentiles_;  
};

template<class TContainer1D>
Real mean( ITContainer1D<Real, TContainer1D> const&  V)
{
  // no samples
  if (V.empty())
  { return Arithmetic<Real>::NA();}

  // get dimensions
  Integer nobs = V.size();
  // sum the samples
  Real sum  = 0.0;
  for (Integer i=V.first(); i<=V.last(); i++)
  { if (Arithmetic<Real>::isFinite(V[i]))
      sum += V[i];
    else
      nobs--;
  }
  // compute the mean
  if (nobs) return sum /= (Real)nobs;
  else      return Arithmetic<Real>::NA();
}

template<class TContainer1D>
Real mean( ITContainer1D<Real, TContainer1D> const&  V
         , ITContainer1D<Real, TContainer1D> const& W
         )
{
  // no samples
  if (V.empty()) { return Arithmetic<Real>::NA();}

  // get dimensions
  Real sum  = 0.0;
  // if the weight are not of the same size, ignore them
  if (V.range() != W.range()) return mean(V);

  // sum the weighted samples
  Real nweight = 0.0;
  for (Integer i=V.first(); i<=V.last(); i++)
  { if ( (Arithmetic<Real>::isFinite(V[i]))
       &&(Arithmetic<Real>::isFinite(W[i]))
       )
    {
      Real weight  = abs(W[i]);
      nweight  += weight;
      sum      += weight * V[i];
    }
  }
  // compute the mean
  return (nweight) ? sum /= nweight : 0.;
}

template<class TContainer1D>
Real varianceWithFixedMean( ITContainer1D<Real, TContainer1D> const& V
                          , Real const& mu
                          , bool unbiased = false
                          )
{
  // no samples
  if (V.empty()) { return Arithmetic<Real>::NA();}

  // get dimensions
  const Integer  first = V.first(), last = V.last();
  Integer nobs = V.size();
  // sum
  Real sum  = 0.0, var  = 0.0, dev;
  for (Integer i=first; i<=last; i++)
  {
    if (Arithmetic<Real>::isFinite(V[i]))
    {
      sum += (dev = V[i] - mu); // deviation from the mean
      var += (dev*dev);         // squared value
    }
    else nobs--;
  }
  // compute the variance
  if (unbiased)
  {
    return (nobs > 1) ? (var - (sum*sum)/(Real)nobs)/(Real)(nobs -1)
                      : Arithmetic<Real>::NA();

  }
  // variance
  return (nobs > 0) ? (var - (sum*sum)/(Real)nobs)/(Real)(nobs)
                    : Arithmetic<Real>::NA();
}

template<class TContainer1D>
Real varianceWithFixedMean( ITContainer1D<Real, TContainer1D> const& V
                          , ITContainer1D<Real, TContainer1D> const& W
                          , Real const& mu
                          , bool unbiased = false
                          )
{
  // no samples
  if (V.empty()) { return Arithmetic<Real>::NA();}

  // if the weight are not of the same size, ignore them
  if (V.range() != W.range()) return varianceWithFixedMean(V, mu);

  // get dimensions
  const Integer  first = V.first(), last = V.last();
  // sum
  Real dev, sum = 0.0, var = 0.0, nweight = 0.0, nweight2 = 0.0;
  for (Integer i=first; i<=last; i++)
  { if ( Arithmetic<Real>::isFinite(V[i]) && Arithmetic<Real>::isFinite(W[i]) )
    {
      Real weight = abs(W[i]);
      nweight    += weight;
      nweight2   += weight * weight;
      sum        += weight*(dev = V[i]-mu); // deviation from the mean
      var        += weight*(dev*dev);       // squared value
    }
  }
  // compute the variance
  if (unbiased)
  {
    return (nweight*nweight - nweight2 > 0.) ? (var - sum*sum/nweight)/(nweight - nweight2/nweight)
                                             : Arithmetic<Real>::NA();

  }
  return (nweight) ? (var - sum*sum)/(nweight) : Arithmetic<Real>::NA();
}

template<class TContainer1D>
Real variance( ITContainer1D<Real, TContainer1D> const& V, bool unbiased = false)
{
    // no samples
    if (V.empty()) { return Arithmetic<Real>::NA();}
    
    // Compute the mean
    Real mu = mean<TContainer1D>(V);
    
    // get dimensions
    const Integer  first = V.first(), last = V.last();
    Integer nobs = V.size();
    // sum
    Real sum  = 0.0, var  = 0.0, dev;
    for (Integer i=first; i<=last; i++)
    {
        if (Arithmetic<Real>::isFinite(V[i]))
        {
            sum += (dev = V[i] - mu); // deviation from the mean
            var += (dev*dev);         // squared value
        }
        else nobs--;
    }
    // compute the variance
    if (unbiased)
    {
        return (nobs > 1) ? (var - (sum*sum)/(Real)nobs)/(Real)(nobs -1)
        : Arithmetic<Real>::NA();
        
    }
    return (nobs > 0) ? (var - (sum*sum)/(Real)nobs)/(Real)(nobs)
    : Arithmetic<Real>::NA();
}

template<class TContainer1D>
Real variance( ITContainer1D<Real, TContainer1D> const& V
              , ITContainer1D<Real, TContainer1D> const& W
              , bool unbiased = false
              )
{
    // no samples
    if (V.empty()) { return Arithmetic<Real>::NA();}
    
    // Compute the mean if necessary
    Real mu = mean(V, W);
    
    // if the weight are not of the same size, ignore them
    if (V.range() != W.range()) return varianceWithFixedMean(V, mu);
    
    // get dimensions
    const Integer  first = V.first(), last = V.last();
    // sum
    Real dev, sum = 0.0, var = 0.0, nweight = 0.0, nweight2 = 0.0;
    for (Integer i=first; i<=last; i++)
    { if ( Arithmetic<Real>::isFinite(V[i]) && Arithmetic<Real>::isFinite(W[i]) )
    {
        Real weight = abs(W[i]);
        nweight    += weight;
        nweight2   += weight * weight;
        sum        += weight*(dev = V[i]-mu); // deviation from the mean
        var        += weight*(dev*dev);       // squared value
    }
    }
    // compute the variance
    if (unbiased)
    {
        return (nweight*nweight - nweight2 > 0.) ? (var - sum*sum/nweight)/(nweight - nweight2/nweight)
        : Arithmetic<Real>::NA();
        
    }
    return (nweight) ? (var - sum*sum)/(nweight) : Arithmetic<Real>::NA();
}
    

}  // namespace Stat

}  // namespace STK

#endif /*STK_STAT_UNIVARIATEREAL_H*/