statistics.hpp

Go to the documentation of this file.

/*
 * Copyright(c):
 * Signal Image and Communications (SIC) Department 
 * http://www.sic.sp2mi.univ-poitiers.fr/ 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - XLIM  Institute UMR CNRS 7252  http://www.xlim.fr/ 
 * 
 * and 
 * 
 * D2 Fluid, Thermic and Combustion 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - PPRIME Institute -  UPR CNRS 3346  http://www.pprime.fr
 *  - ISAE-ENSMA http://www.ensma.fr
 *
 * Contributor(s):
 * The SLIP team,
 * Benoit Tremblais <tremblais_AT_sic.univ-poitiers.fr>,
 * Laurent David <laurent.david_AT_lea.univ-poitiers.fr>, 
 * Ludovic Chatellier <ludovic.chatellier_AT_univ-poitiers.fr>, 
 * Lionel Thomas <lionel.thomas_AT_univ-poitiers.fr>, 
 * Denis Arrivault <arrivault_AT_sic.univ-poitiers.fr>, 
 * Julien Dombre <julien.dombre_AT_univ-poitiers.fr>.
 *
 * Description:
 * The Simple Library of Image Processing (SLIP) is a new image processing 
 * library. It is written in the C++ language following as much as possible 
 * the ISO/ANSI C++ standard. It is consequently compatible with any system  
 * satisfying the ANSI C++ complience. It works on different Unix , Linux ,  
 * Mircrosoft Windows and Mac OS X plateforms. SLIP is a research library that  
 * was created by the Signal, Image and Communications (SIC) departement of  
 * the XLIM, UMR 7252 CNRS Institute in collaboration with the Fluids, Thermic  
 * and Combustion departement of the P', UPR 3346 CNRS Institute of the  
 * University of Poitiers.
 * 
 * The SLIP Library source code has been registered to the APP (French Agency 
 * for the Protection of Programs) by the University of Poitiers and CNRS, 
 * under  registration number IDDN.FR.001.300034.000.S.P.2010.000.21000.

 * http://www.sic.sp2mi.univ-poitiers.fr/slip/
 *
 * This software is governed by the CeCILL-C license under French law and
 * abiding by the rules of distribution of free software.  You can  use, 
 * modify and/ or redistribute the software under the terms of the CeCILL-C 
 * license as circulated by CEA, CNRS and INRIA at the following URL 
 * http://www.cecill.info.
 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.
 * 
 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the 
 * same conditions as regards security. 
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL-C license and that you accept its terms.
 */



#ifndef SLIP_STATISTICS_HPP
#define SLIP_STATISTICS_HPP

#include <iostream>
#include <algorithm>
#include <iterator>
#include <numeric>
#include <cassert>
#include <cmath>
#include <set>
#include <vector>
#include "Array.hpp"
#include "macros.hpp"
#include "iterator_types.hpp"
#include "Block.hpp"

namespace slip
{


template <typename T>
struct Statistics
{
  T& min() 
  { 
    return stat_[0];
  }
 
  T& first_quartile()
  {
    return stat_[1];
  }

  T& median()
  {
    return stat_[2];
  }

  T& third_quartile()
  {
    return stat_[3];
  }

  T& max() 
  { 
    return stat_[4];
  }

  T& mean()
  { 
    return stat_[5];
  }

  T& std_dev()
  { 
    return stat_[6];
  }

  T& skewness()
  { 
    return stat_[7];
  }

  T& kurtosis()
  {
    return stat_[8];
  }

  T& cardinal()
  {
    return stat_[9];
  }
  
  slip::block<T,10> all()
  {
    return stat_;
  }

  slip::block<T,10> stat_;
};  

  template<typename Integer, typename>
  struct __cardinal
  {
    template <typename _II>
    static Integer
    cardinal(_II first, _II last)
    {
      return static_cast<Integer>(last-first);
    
    }
  };

  template<typename Integer>
  struct __cardinal<Integer,std::random_access_iterator_tag>
  {
    template <typename _II>
    static Integer
    cardinal(_II first, _II last)
    {
      return static_cast<Integer>(last-first);
    }
  };

  template<typename Integer>
  struct __cardinal<Integer,std::random_access_iterator2d_tag>
  {
    template <typename _II>
    static Integer
    cardinal(_II first, _II last)
    {
      Integer count = 0;
      while(first != last)
      {
                count++;
                first++;
      }
      return count;

    }
  };


   template<typename Integer>
  struct __cardinal<Integer,std::random_access_iterator3d_tag>
  {
    template <typename _II>
    static Integer
    cardinal(_II first, _II last)
    {
      return __cardinal<Integer,std::random_access_iterator2d_tag>::cardinal(first,last);
    }
  };

  template<typename Integer,typename InputIterator>
  inline
  Integer cardinal(InputIterator first, 
                                   InputIterator last)
  {
    typedef typename std::iterator_traits<InputIterator>::iterator_category _Category;
  
    return __cardinal<Integer,_Category>::cardinal(first,last);
  }

  template<typename Value_T,typename InputIterator>
  inline
  Value_T mean(InputIterator first, 
                       InputIterator last)
  {
    assert(first != last);
    Value_T sum = Value_T(0);
    long int count = 0;
    while(first != last)
      {
                sum += *first++;
                ++count;
      }
    return Value_T(sum) / Value_T(count);
  }

  template<typename Value_T,typename InputIterator>
  inline
  Value_T weighted_mean(InputIterator first, 
                                                InputIterator last,
                                                InputIterator w_first,
                                                InputIterator w_last)
  {
    assert(first != last);
    assert(slip::cardinal<std::size_t>(first,last) == 
                   slip::cardinal<std::size_t>(w_first,w_last));

    Value_T proj = Value_T(0);
    Value_T sum  = Value_T(0);
    for (; first != last; ++first, ++w_first)
       {
                 proj = proj + (*first) * (*w_first);
                 sum  = sum  + *w_first;
       }
    return Value_T(proj) / Value_T(sum);
  }

 

 

  
  template<typename T, int N, typename InputIterator>
  inline
  T nth_moment(InputIterator first, InputIterator last, T mean)
  {
    assert(last != first);
    T __init = T();
    int count = 0;
     for (; first != last; ++first)
       {
                __init = __init + slip::nth_power<N>(*first - mean);
                count++;
       }
     return __init / count;
  }

 
  template<typename T, typename InputIterator, typename InputIterator2>
  inline
  T covariance(InputIterator first, InputIterator last, 
                       InputIterator2 first2, 
                       T mean1, T mean2)
  {
    assert(last != first);
    T __init = T();
    long int count = 0;
     for (; first != last; ++first, ++first2)
       {
                __init = __init + (*first - mean1) * (*first2 - mean2);
                count++;
       }
     return __init / count;
  }

  template<typename T, typename InputIterator, typename InputIterator2>
  inline
  T unbiased_covariance(InputIterator first, InputIterator last, 
                                                InputIterator2 first2, 
                                                T mean1, T mean2)
  {
    assert(last != first);
    T __init = T();
    long int count = 0;
     for (; first != last; ++first, ++first2)
       {
                __init = __init + (*first - mean1) * (*first2 - mean2);
                count++;
       }
     T result = T(); 
     if(count != 1)
       {
                 result = __init / (count - 1);
       }
   
     return result;
  }
  template<typename T, typename InputIterator>
  inline
  T variance(InputIterator first, InputIterator last, T mean)
  {
    return slip::covariance(first,last,first,mean,mean);
  }

  template<typename T, typename InputIterator>
  inline
  T unbiased_variance(InputIterator first, InputIterator last, T mean)
  {
    return slip::unbiased_covariance(first,last,first,mean,mean);
  }

  template<typename T, typename InputIterator>
  inline
  T std_dev(InputIterator first, InputIterator last, T mean)
  {
    return std::sqrt(slip::variance(first,last,mean));
  }

  template<typename T, typename InputIterator>
  inline
  T unbiased_std_dev(InputIterator first, InputIterator last, T mean)
  {
    return std::sqrt(slip::unbiased_variance(first,last,mean));
  }

 
  template<typename T, typename InputIterator>
  inline
  T kurtosis(InputIterator first, InputIterator last, T mean)
  {
    assert(first != last);
    T __init  = T(0);
    T __init2 = T(0);
    long int count = 0;
     for (; first != last; ++first)
       {
                 T p2 = slip::nth_power<2>(*first - mean);
                 __init  = __init  + p2;
                 __init2 = __init2 + (p2 * p2);
                 count++;
       }
     T result = T(0);
     if(__init != T(0))
       {
                 result = (__init2 * T(count)) / (__init * __init);
       }
     return result;
  }
 

  template<typename T, typename InputIterator>
  inline
  T unbiased_kurtosis(InputIterator first, InputIterator last, T mean)
  {
    assert(slip::cardinal<int>(first,last) > 3);
    T __init  = T(0);
    T __init2 = T(0);
    long int count = 0;
    for (; first != last; ++first)
       {
                 T p2 = slip::nth_power<2>(*first - mean);
                 __init  = __init  + p2;
                 __init2 = __init2 + (p2 * p2);
                 count++;
       }
    T count_m1 = T(count - 1);
    T count_m2 = T(count - 2);
    T count_m3 = T(count - 3);
  
    T result = T(0);
    if(__init != T(0))
      {
    
                result = (((T(count+1)*T(count)* count_m1) * __init2 ) 
                                  / ( (count_m2 * count_m3) * (__init * __init) ))
                  - ( (T(3)* (count_m1*count_m1)) / (count_m2 * count_m3)) + T(3);
      }
  return result;
  }

  template<typename T, typename InputIterator>
  inline
  T skewness(InputIterator first, InputIterator last, T mean)
  {
    assert(first != last);
    T __init  = T();
    T __init2 = T();
    long int count = 0;
     for (; first != last; ++first)
       {
                 T p1 = (*first - mean);
                 T p2 =  p1 * p1;
                 __init  = __init  + p2;
                 __init2 = __init2 + (p2 * p1);
                 count++;
       }
     T result = T(0);
     if(__init != T(0))
       {
                 result = __init2 / (T(count)*slip::nth_power<3>(std::sqrt(__init/T(count))));
       }
     return result;
  }

  template<typename T, typename InputIterator>
  inline
  T unbiased_skewness(InputIterator first, InputIterator last, T mean)
  {
    assert(slip::cardinal<int>(first,last) > 2);
    T __init  = T();
    T __init2 = T();
    long int count = 0;
     for (; first != last; ++first)
       {
                 T p1 = (*first - mean);
                 T p2 =  p1 * p1;
                 __init  = __init  + p2;
                 __init2 = __init2 + (p2 * p1);
                 count++;
       }
    T result = T(0);
    if(__init != T(0))
      {
                result = (__init2 * T(count)) / (T(count-1)*T(count-2)*slip::nth_power<3>(std::sqrt(__init/T(count-1))));
      }
    return result;
}
 
 
  template<typename T, typename InputIterator>
  inline
  T rms(InputIterator first, InputIterator last)
  {
    return std::sqrt(slip::covariance(first,last,first,T(0),T(0)));
  }


  template<typename InputIterator>
  inline
  typename std::iterator_traits<InputIterator>::value_type
  median_from_sorted_data(InputIterator first, InputIterator last)
  {
    return *(first + (last - first) / 2); 
  }
  

  template<typename RandomAccessIterator, typename Size>
  inline
  typename std::iterator_traits<RandomAccessIterator>::value_type
  median_from_sorted_data_n(RandomAccessIterator first, Size n)
  {
    return *(first + n/2); 
  }

  

 

  template<typename RandomAccessIterator, typename Size>
  inline
  typename std::iterator_traits<RandomAccessIterator>::value_type
  median_from_data_n(RandomAccessIterator first, RandomAccessIterator last, 
                                     Size n)
  {
    assert(last != first);
    slip::Array<typename std::iterator_traits<RandomAccessIterator>::value_type> 
      tmp(n,first,last);
    std::sort(tmp.begin(),tmp.end());
    return slip::median_from_sorted_data_n(tmp.begin(),n); 
  }

  template<typename RandomAccessIterator, 
                   typename Size, 
                   typename StrictWeakOrdering>
  inline
  typename std::iterator_traits<RandomAccessIterator>::value_type
  median_from_data_n(RandomAccessIterator first, RandomAccessIterator last, 
                                     Size n, 
                                     StrictWeakOrdering comp)
  {
    assert(last != first);
    slip::Array<typename std::iterator_traits<RandomAccessIterator>::value_type> 
      tmp(n,first,last);
    std::sort(tmp.begin(),tmp.end(),comp);
    return slip::median_from_sorted_data_n(tmp.begin(),n); 
  }


  template<typename>
  struct __median
  {
    template <typename _II>
    static typename std::iterator_traits<_II>::value_type
    median(_II first, _II last)
    {
      typedef typename std::iterator_traits<_II>::value_type value_type;
      return slip::median_from_data_n(first,last,slip::cardinal<int>(first,last));
    }
    template <typename _II, typename StrictWeakOrdering>
    static typename std::iterator_traits<_II>::value_type
    median(_II first, _II last, StrictWeakOrdering comp)
    {
      typedef typename std::iterator_traits<_II>::value_type value_type;
      return slip::median_from_data_n(first,last,slip::cardinal<int>(first,last),comp);
    }
    template <typename _II,typename SizeType>
    static typename std::iterator_traits<_II>::value_type
    median_n(_II first, _II last, SizeType n)
    {
      typedef typename std::iterator_traits<_II>::value_type value_type;
      return slip::median_from_data_n(first,last,n);
    }
    template <typename _II,typename SizeType, typename StrictWeakOrdering>
    static typename std::iterator_traits<_II>::value_type
    median_n(_II first, _II last, SizeType n, StrictWeakOrdering comp)
    {
      typedef typename std::iterator_traits<_II>::value_type value_type;
      return slip::median_from_data_n(first,last,n,comp);
    }
  };

 
  template<>
  struct __median<std::random_access_iterator_tag>
  {
    template <typename _II>
    static typename std::iterator_traits<_II>::value_type
    median(_II first, _II last)
    {
      typedef typename std::iterator_traits<_II>::value_type value_type;
      return slip::median_from_data_n(first,last,(last-first));
    }
    template <typename _II, typename StrictWeakOrdering>
    static typename std::iterator_traits<_II>::value_type
    median(_II first, _II last, StrictWeakOrdering comp)
    {
      typedef typename std::iterator_traits<_II>::value_type value_type;
      return slip::median_from_data_n(first,last,(last-first),comp);
    }
     template <typename _II,typename SizeType>
    static typename std::iterator_traits<_II>::value_type
    median_n(_II first, _II last, SizeType n)
    {
      typedef typename std::iterator_traits<_II>::value_type value_type;
      return slip::median_from_data_n(first,last,n);
    }

    template <typename _II,typename SizeType, typename StrictWeakOrdering>
    static typename std::iterator_traits<_II>::value_type
    median_n(_II first, _II last, SizeType n, StrictWeakOrdering comp)
    {
      typedef typename std::iterator_traits<_II>::value_type value_type;
      return slip::median_from_data_n(first,last,n,comp);
    }
  };
  
  template<typename RandomAccessIterator>
  inline
  typename std::iterator_traits<RandomAccessIterator>::value_type
  median(RandomAccessIterator first, RandomAccessIterator last)
  {
     typedef typename std::iterator_traits<RandomAccessIterator>::iterator_category _Category;
  
     return __median<_Category>::median(first,last);
  }


  template<typename RandomAccessIterator, typename StrictWeakOrdering>
  inline
  typename std::iterator_traits<RandomAccessIterator>::value_type
  median(RandomAccessIterator first, RandomAccessIterator last, StrictWeakOrdering comp)
  {
     typedef typename std::iterator_traits<RandomAccessIterator>::iterator_category _Category;
  
     return __median<_Category>::median(first,last,comp);
  }


  template<typename RandomAccessIterator, typename SizeType>
  inline
  typename std::iterator_traits<RandomAccessIterator>::value_type
  median_n(RandomAccessIterator first, RandomAccessIterator last, SizeType n)
  {
     typedef typename std::iterator_traits<RandomAccessIterator>::iterator_category _Category;
  
     return __median<_Category>::median_n(first,last,n);
  }
  
  template<typename RandomAccessIterator, typename SizeType, typename StrictWeakOrdering>
  inline
  typename std::iterator_traits<RandomAccessIterator>::value_type
  median_n(RandomAccessIterator first, RandomAccessIterator last, SizeType n, StrictWeakOrdering comp)
  {
     typedef typename std::iterator_traits<RandomAccessIterator>::iterator_category _Category;
  
     return __median<_Category>::median_n(first,last,n,comp);
  }
  

  template<typename RandomAccessIterator>
  inline
  typename std::iterator_traits<RandomAccessIterator>::value_type
  first_quartile(RandomAccessIterator first, RandomAccessIterator last)
  {
    assert(last != first);
    std::size_t n = slip::cardinal<std::size_t>(first,last);
    slip::Array<typename std::iterator_traits<RandomAccessIterator>::value_type> 
      tmp(n,first,last);
    std::sort(tmp.begin(),tmp.end());
    return *(tmp.begin()+n/4);
  }

  template<typename RandomAccessIterator>
  inline
  typename std::iterator_traits<RandomAccessIterator>::value_type
  third_quartile(RandomAccessIterator first, RandomAccessIterator last)
  {
    assert(last != first);
    std::size_t n=slip::cardinal<std::size_t>(first,last);
    slip::Array<typename std::iterator_traits<RandomAccessIterator>::value_type> 
      tmp(n,first,last);
    std::sort(tmp.begin(),tmp.end());
    return *((tmp.end()-1)-n/4);
  }
  


  template<typename ForwardIterator>
  void n_max_elements(ForwardIterator first, ForwardIterator last,
                                      std::vector<ForwardIterator>& max, 
                                      const std::size_t n = 1)
  {
    
    assert(n <= slip::cardinal<std::size_t>(first,last));
    assert(max.size() == n);
    assert(n > 0);
    if(n == 1)
      {
                max[0] = std::max_element(first,last);
      }
    else
      {
                std::multiset<ForwardIterator,slip::gt_it<ForwardIterator> > S;
                
                while(first != last)
                  {
                    S.insert(first);
                    first++;
                  }
                
                typename std::multiset<ForwardIterator,gt_it<ForwardIterator> >::iterator it = S.begin();
                for(std::size_t i = 0; i < n; ++i)
                  {
                    max[i] = *it++;
                  }
                
      }
}


  template<typename ForwardIterator>
  void n_min_elements(ForwardIterator first, ForwardIterator last,
                                      std::vector<ForwardIterator>& min, 
                                      const std::size_t n = 1)
  {
    assert(n <= slip::cardinal<std::size_t>(first,last));
    assert(min.size() == n);
    assert(n > 0);
    if(n == 1)
      {
                min[0] = std::min_element(first,last);
      }
    else
      {
                std::multiset<ForwardIterator,slip::lt_it<ForwardIterator> > S;
                
                while(first != last)
                  {
                    S.insert(first);
                    first++;
                  }
                
                typename std::multiset<ForwardIterator,lt_it<ForwardIterator> >::iterator it = S.begin();
                for(std::size_t i = 0; i < n; ++i)
                  {
                    min[i] = *it++;
                  }
                
      }
  }

  template<typename RandomAccessIterator, typename RandomAccessIterator2>
  void center(RandomAccessIterator first, RandomAccessIterator last,
                      RandomAccessIterator2 out_first)
  { 
    typedef typename std::iterator_traits<RandomAccessIterator2>::value_type value_type;
    value_type mean  = slip::mean<value_type>(first,last);
    
    for (; first != last; ++first, ++out_first)
      {
                *out_first = (*first - mean);
      }
  }

  template<typename RandomAccessIterator, typename RandomAccessIterator2>
  void studentize(RandomAccessIterator first, RandomAccessIterator last,
                                  RandomAccessIterator2 out_first)
  { 
    typedef typename std::iterator_traits<RandomAccessIterator2>::value_type value_type;
    value_type mean  = slip::mean<value_type>(first,last);
    value_type std_dev = slip::std_dev<value_type>(first,last,mean);
    if(std_dev != value_type(0))
      {
                for (; first != last; ++first, ++out_first)
                  {
                    *out_first = (*first - mean) / std_dev;
                  }
      }
    else
      {
                for (; first != last; ++first, ++out_first)
                  {
                    *out_first = (*first - mean);
                  }
      }
  }
  



   template<typename Value_T,
                    typename InputIterator, 
                    typename MaskIterator>
   inline
   Value_T mean_mask(InputIterator first, InputIterator last, 
                                     MaskIterator mask_first, 
                                     typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
   {
    assert(first != last); 
    Value_T sum = Value_T(0); 
    long int count = 0; 
    while(first != last)
        {
            if(*mask_first == value)
            { 
                                sum += *first; 
                                ++count; 
                    }

        mask_first++; 
        first++;
        }
    Value_T result = Value_T(0);
    if (count != 0)
      {
                result = Value_T(sum) / Value_T(count);
      }
    return result;
   }

  template<typename T, int N, typename InputIterator, typename MaskIterator>
  inline
  T nth_moment_mask(InputIterator first, InputIterator last, MaskIterator mask_first, T mean, typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    T __init = T();
    int count = 0;
     for (; first != last; ++first,++mask_first)
       {
                if(*mask_first == value)
            { 
                                __init = __init + slip::nth_power<N>(*first - mean);
                                count++;
                     }
       }
     T result = T(0);
     if(count != 0)
       {
                 result = __init / static_cast<T>(count);
       }
     return result;
  }

 
  template<typename T, typename InputIterator, typename InputIterator2, typename MaskIterator>
  inline
  T covariance_mask(InputIterator first, InputIterator last, 
                       InputIterator2 first2, MaskIterator mask_first, 
                                    T mean1, T mean2, typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    assert(first != last);
    T __init = T();
    long int count = 0;
     for (; first != last; ++first, ++first2,++mask_first)
       {
                if(*mask_first == value)
            { 
                                __init = __init + (*first - mean1) * (*first2 - mean2);
                                count++;
                    }
       }
     T result = T(0);
     if(count != 0)
       {
                 result = __init / static_cast<T>(count);
       }
     return result;
  }

  template<typename T, typename InputIterator, typename InputIterator2, typename MaskIterator>
  inline
  T unbiased_covariance_mask(InputIterator first, InputIterator last, 
                       InputIterator2 first2, MaskIterator mask_first, 
                                    T mean1, T mean2, typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    T __init = T();
    long int count = 0;
     for (; first != last; ++first, ++first2,++mask_first)
       {
                if(*mask_first == value)
            { 
                                __init = __init + (*first - mean1) * (*first2 - mean2);
                                count++;
                    }
       }

     T result = T(0);
     if(count != 0)
       {
                 result = __init / static_cast<T>(count-1);
       }
     return result;
  }


  template<typename T, typename InputIterator, typename MaskIterator>
  inline
  T variance_mask(InputIterator first, InputIterator last, MaskIterator mask_first, T mean, typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    return slip::covariance_mask(first,last,first,mask_first, mean,mean,value);
  }


  template<typename T, typename InputIterator, typename MaskIterator>
  inline
  T unbiased_variance_mask(InputIterator first, InputIterator last, MaskIterator mask_first, T mean, typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    return slip::unbiased_covariance_mask(first,last,first,mask_first, mean,mean,value);
  }

  template<typename T, typename InputIterator, typename MaskIterator>
  inline
  T std_dev_mask(InputIterator first, InputIterator last, MaskIterator mask_first, T mean, typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    return std::sqrt(slip::variance_mask(first,last,mask_first,mean,value));
  }


  template<typename T, typename InputIterator, typename MaskIterator>
  inline
  T unbiased_std_dev_mask(InputIterator first, InputIterator last, MaskIterator mask_first, T mean, typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    return std::sqrt(slip::unbiased_variance_mask(first,last,mask_first,mean,value));
  }

 
  template<typename T, typename InputIterator, typename MaskIterator>
  inline
  T kurtosis_mask(InputIterator first, InputIterator last, MaskIterator mask_first, T mean, typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    T __init  = T(0);
    T __init2 = T(0);
    long int count = 0;
     for (; first != last; ++first,++mask_first)
       {
                if (*mask_first == value)
                  {
                    T p2 = slip::nth_power<2>(*first - mean);
                    __init  = __init  + p2;
                    __init2 = __init2 + (p2 * p2);
                    count++;
                  }
       }
     T result = T(0);
     if( __init != T(0))
       {
                 result =  (__init2 * count) / (__init * __init);
       }
     return result;
  }
 

  template<typename T, 
                   typename InputIterator, 
                   typename MaskIterator>
  inline
  T unbiased_kurtosis_mask(InputIterator first, InputIterator last, 
                                                   MaskIterator mask_first, 
                                                   T mean, 
                                                   typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    T __init  = T(0);
    T __init2 = T(0);
    long int count = 0;
     for (; first != last; ++first,++mask_first)
       {
                if (*mask_first == value)
                  {
                    T p2 = slip::nth_power<2>(*first - mean);
                    __init  = __init  + p2;
                    __init2 = __init2 + (p2 * p2);
                    count++;
                  }
       }

     T count_m1 = T(count - 1);
     T count_m2 = T(count - 2);
     T count_m3 = T(count - 3);
     assert(count > 3);
     T result = T(0);
     if(__init != T(0))
       {
                 result = (((T(count+1)*T(count)* count_m1) * __init2 ) 
                                   / ( (count_m2 * count_m3) * (__init * __init) ))
                   - ( (T(3)* (count_m1*count_m1)) / (count_m2 * count_m3)) + T(3);
       }
     return result;
     
  }
 

  template<typename T, 
                   typename InputIterator, 
                   typename MaskIterator>
  inline
  T skewness_mask(InputIterator first, InputIterator last, 
                                  MaskIterator mask_first, 
                                  T mean, 
                                  typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    T __init  = T(0);
    T __init2 = T(0);
    long int count = 0;
     for (; first != last; ++first,++mask_first)
       {
                if (*mask_first == value)
                  {
                    T p1 = (*first - mean);
                    T p2 =  p1 * p1;
                    __init  = __init  + p2;
                    __init2 = __init2 + (p2 * p1);
                    count++;
                 
                  }
       }
     assert(count > 0);
     T result = T(0);
     if(__init != T(0))
       {
                 result = __init2 / (T(count)*slip::nth_power<3>(std::sqrt(__init/T(count))));
       }
     return result;
  }
 

  template<typename T, 
                   typename InputIterator, 
                   typename MaskIterator>
  inline
  T unbiased_skewness_mask(InputIterator first, InputIterator last, 
                                                   MaskIterator mask_first, 
                                                   T mean, 
                                                   typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    T __init  = T(0);
    T __init2 = T(0);
    long int count = 0;
     for (; first != last; ++first,++mask_first)
       {
                if (*mask_first == value)
                  {
                    T p1 = (*first - mean);
                    T p2 =  p1 * p1;
                    __init  = __init  + p2;
                    __init2 = __init2 + (p2 * p1);
                    count++;
                 
                  }
       }
     assert(count > 2);
     T result = T(0);
     if(__init != T(0))
       {
                 result = (__init2 * T(count)) / (T(count-1)*T(count-2)*slip::nth_power<3>(std::sqrt(__init/T(count-1))));
       }
     return result;
}

  template<typename T, 
                   typename InputIterator, 
                   typename MaskIterator>
  inline
  T rms_mask(InputIterator first, InputIterator last, 
                     MaskIterator mask_first,
                     typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    return std::sqrt(slip::covariance_mask(first,last,first,mask_first,         T(0),T(0),value));
  }

  template<typename Integer,typename MaskIterator>
  inline
  Integer cardinal_mask(MaskIterator mask_first, 
                                                MaskIterator mask_last,
                                                typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    Integer count = 0;
    while(mask_first != mask_last)
      {
                if (*mask_first == value)
                  {
                    count++;
                  }
                mask_first++;
      }
    return count;
  }

  template<typename InputIterator, 
                   typename MaskIterator>
  inline
  typename std::iterator_traits<InputIterator>::value_type
  median_from_sorted_data_mask(InputIterator first, InputIterator last, 
                                                       MaskIterator mask_first, 
                                                       typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    MaskIterator mask_last=mask_first+(last - first);
    std::size_t mask_count=slip::cardinal_mask<std::size_t>(mask_first,mask_last,value);
    assert(mask_count != 0);
    std::size_t count = 0;
    while(first != last)
      {
                if (*mask_first == value)
                  {
                                count++;
                  }
                if (count>mask_count/2)
                  {
                                return *first;
                  }
                first++;
        mask_first++;
      }    
    std::cerr << "Error computing median_from_sorted_data_mask"<<std::endl;
    return *first; 
  }
  
  template<typename InputIterator, typename MaskIterator>
  inline
  typename std::iterator_traits<InputIterator>::value_type
  median_mask(InputIterator first, InputIterator last, 
                      MaskIterator mask_first, 
                      typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    std::size_t mask_count = 
      slip::cardinal_mask<std::size_t>(mask_first,mask_first+(last-first),
                                                                       value);

    assert(mask_count != 0);
    slip::Array<typename std::iterator_traits<InputIterator>::value_type> tmp(mask_count); 
    std::size_t count = 0;
    while(first != last)
      {
                if (*mask_first == value)
                  {
                    tmp[count]=*first;
                    count++;
                  }
                first++;
        mask_first++;
      }    
    
    std::sort(tmp.begin(),tmp.end());
    return slip::median_from_sorted_data_n(tmp.begin(),mask_count); 

  }

  template<typename InputIterator, 
                   typename MaskIterator>
  inline
  typename std::iterator_traits<InputIterator>::value_type
  first_quartile_mask(InputIterator first, InputIterator last, 
                                      MaskIterator mask_first, 
                                      typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    std::size_t mask_count = slip::cardinal_mask<std::size_t>(mask_first,mask_first+(last-first),value);
    assert(mask_count != 0);
    slip::Array<typename std::iterator_traits<InputIterator>::value_type> tmp(mask_count);
    std::size_t count = 0;
    while(first != last)
      {
                if (*mask_first == value)
                  {
                    tmp[count]=*first;
                    count++;
                  }
                first++;
        mask_first++;

      }    

    std::sort(tmp.begin(),tmp.end());
    return slip::first_quartile(tmp.begin(),tmp.end()); 

  }

  template<typename InputIterator, typename MaskIterator>
  inline
  typename std::iterator_traits<InputIterator>::value_type
  third_quartile_mask(InputIterator first, InputIterator last, MaskIterator mask_first, typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
  {
    std::size_t mask_count=slip::cardinal_mask<std::size_t>(mask_first,mask_first+(last-first),value);
    assert(mask_count != 0);
    slip::Array<typename std::iterator_traits<InputIterator>::value_type> tmp(mask_count);

    std::size_t count = 0;
    while(first != last)
      {
                if (*mask_first == value)
                  {
                    tmp[count]=*first;
                    count++;
                  }
                first++;
        mask_first++;
      }    

    std::sort(tmp.begin(),tmp.end());
    return slip::third_quartile(tmp.begin(),tmp.end()); 

  }

  template<typename InputIterator, typename MaskIterator, typename OutputIterator>
   inline
   void center_mask(InputIterator first, InputIterator last, 
                                    MaskIterator mask_first, 
                                    OutputIterator out_first, 
                                    typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
   {
    assert(first != last); 
    typedef typename std::iterator_traits<OutputIterator>::value_type value_type;
 
    value_type mean = slip::mean_mask<value_type>(first,last,mask_first,value);
    
    for (; first != last; ++first, ++out_first, ++mask_first)
      { 
                if(*mask_first == value)
                  {
                    *out_first = (*first - mean);
                  }
      }
                 
   }


  template<typename InputIterator, typename MaskIterator, typename OutputIterator>
   inline
   void studentize_mask(InputIterator first, InputIterator last, 
                                    MaskIterator mask_first, 
                                    OutputIterator out_first, 
                                    typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
   {
    assert(first != last); 
    typedef typename std::iterator_traits<OutputIterator>::value_type value_type;
 
    value_type mean = slip::mean_mask<value_type>(first,last,mask_first,value);
    value_type std_dev = slip::std_dev_mask<value_type>(first,last,mask_first,mean,value);
    
    if(std_dev != value_type(0))
      {
                for (; first != last; ++first, ++out_first, ++mask_first)
                  { 
                    if(*mask_first == value)
                      {
                                *out_first = (*first - mean) / std_dev;
                      }
                  }
      }
    else
      {
                for (; first != last; ++first, ++out_first, ++mask_first)
                  { 
                    if(*mask_first == value)
                      {
                                *out_first = (*first - mean);
                      }
                  }
      }
                 
   }



  template<typename Integer,typename InputIterator, typename Predicate>
  inline
  Integer cardinal_if(InputIterator first, 
                                   InputIterator last,
                                Predicate pred)
  {
    Integer count = 0;
    while(first != last)
      {
                if (pred(*first))
                  {
                    count++;
                  }
                first++;
      }
    return count;
  }

   template<typename Value_T,typename InputIterator, typename Predicate>
   inline
   Value_T mean_if(InputIterator first, InputIterator last, Predicate pred)
   {
    assert(first != last); 
    Value_T sum = Value_T(0); 
    long int count = 0; 
    while(first != last)
        {
            if(pred(*first))
            { 
                                sum += *first; ++count; 
                    }

        first++;
        }
    Value_T result = Value_T(0);
    if(count != 0)
      {
                result =  Value_T(sum) / Value_T(count);
      }
    return result;
   }

  template<typename T, int N, typename InputIterator, typename Predicate>
  inline
  T nth_moment_if(InputIterator first, InputIterator last, T mean, Predicate pred)
  {
    T __init = T();
    long int count = 0;
     for (; first != last; ++first)
       {
                if(pred(*first))
            { 
                                __init = __init + slip::nth_power<N>(*first - mean);
                                count++;
                     }
       }
     T result = T(0);
     if(count != 0)
       {
                 result = __init / count;
       }
     return result;
  }

 
  template<typename T, 
                   typename InputIterator, 
                   typename InputIterator2, 
                   typename Predicate>
  inline
  T covariance_if(InputIterator first, InputIterator last, 
                                  InputIterator2 first2, 
                                  T mean1, T mean2, Predicate pred)
  {
    T __init = T();
    long int count = 0;
     for (; first != last; ++first, ++first2)
       {
                if(pred(*first))
            { 
                                __init = __init + (*first - mean1) * (*first2 - mean2);
                                count++;
                    }
       }
     T result = T(0);
     if(count != 0)
       {
                 result = __init / count;
       }
     return result;
  }

  template<typename T, typename InputIterator, typename InputIterator2, typename Predicate>
  inline
  T unbiased_covariance_if(InputIterator first, InputIterator last, 
                                  InputIterator2 first2, 
                                  T mean1, T mean2,Predicate pred)
  {
    T __init = T();
    long int count = 0;
     for (; first != last; ++first, ++first2)
       {
                if(pred(*first))
            { 
                                __init = __init + (*first - mean1) * (*first2 - mean2);
                                count++;
                    }
       }
     T result = T(0);
     if(count != 0)
       {
                 result = __init / static_cast<T>(count - 1);
       }
     return result;
  }


  template<typename T, typename InputIterator, typename Predicate>
  inline
  T variance_if(InputIterator first, InputIterator last, T mean, Predicate pred)
  {
    return slip::covariance_if(first,last,first, mean,mean,pred);
  }


  template<typename T, typename InputIterator, typename Predicate>
  inline
  T unbiased_variance_if(InputIterator first, InputIterator last, T mean, Predicate pred)
  {
    return slip::unbiased_covariance_if(first,last,first, mean,mean,pred);
  }

  template<typename T, typename InputIterator, typename Predicate>
  inline
  T std_dev_if(InputIterator first, InputIterator last, T mean, Predicate pred)
  {
    return std::sqrt(slip::variance_if(first,last,mean,pred));
  }


  template<typename T, typename InputIterator, typename Predicate>
  inline
  T unbiased_std_dev_if(InputIterator first, InputIterator last, T mean, Predicate pred)
  {
    return std::sqrt(slip::unbiased_variance_if(first,last,mean,pred));
  }

 
  template<typename T, typename InputIterator, typename Predicate>
  inline
  T kurtosis_if(InputIterator first, InputIterator last, T mean, Predicate pred)
  {
    T __init  = T();
    T __init2 = T();
    long int count = 0;
     for (; first != last; ++first)
       {
                if (pred(*first))
                  {
                                T p2 = slip::nth_power<2>(*first - mean);
                                __init  = __init  + p2;
                                __init2 = __init2 + (p2 * p2);
                                count++;
                  }
       }
     T result = T(0);
     if( __init != T(0))
       {
                 result =  (__init2 * T(count)) / (__init * __init);
       }
     return result;
  }
 

  template<typename T, typename InputIterator, typename Predicate>
  inline
  T unbiased_kurtosis_if(InputIterator first, InputIterator last, T mean, Predicate pred)
  {
    T __init  = T();
    T __init2 = T();
    long int count = 0;
     for (; first != last; ++first)
       {
                if (pred(*first)){
                                T p2 = slip::nth_power<2>(*first - mean);
                                __init  = __init  + p2;
                                __init2 = __init2 + (p2 * p2);
                                count++;
                }
       }
       assert(count > 3);
       T count_m1 = T(count - 1);
       T count_m2 = T(count - 2);
       T count_m3 = T(count - 3);

       T result = T(0);
       if( __init != T(0))
                 {
                   result =   (((T(count+1)*T(count)* count_m1) * __init2 ) 
                                / ( (count_m2 * count_m3) * (__init * __init) ))
                 - ( (T(3)* (count_m1*count_m1)) / (count_m2 * count_m3)) + T(3);
                 }
       return result;
  }


  template<typename T, typename InputIterator, typename Predicate>
  inline
  T skewness_if(InputIterator first, InputIterator last, T mean, Predicate pred)
  {
    T __init  = T();
    T __init2 = T();
    long int count = 0;
     for (; first != last; ++first)
       {
                if (pred(*first))
                  {
                    T p1 = (*first - mean);
                    T p2 =  p1 * p1;
                    __init  = __init  + p2;
                    __init2 = __init2 + (p2 * p1);
                    count++;
                  }
       }
     assert(count > 0);
     T result = T(0);
     if( __init != T(0))
       {
                 result = __init2 / (T(count)*slip::nth_power<3>(std::sqrt(__init/T(count))));
       }
     return result;
  }

  template<typename T, typename InputIterator, typename Predicate>
  inline
  T unbiased_skewness_if(InputIterator first, InputIterator last, T mean, Predicate pred)
  {
    T __init  = T();
    T __init2 = T();
    long int count = 0;
     for (; first != last; ++first)
       {
                if (pred(*first))
                  {
                    T p1 = (*first - mean);
                    T p2 =  p1 * p1;
                    __init  = __init  + p2;
                    __init2 = __init2 + (p2 * p1);
                    count++;
                  }
       }
     assert(count > 2);
     T result = T(0);
     if( __init != T(0))
       {
                 result = (__init2 * T(count)) / (T(count-1)*T(count-2)*slip::nth_power<3>(std::sqrt(__init/T(count-1))));
       }
     return result;
  }

  template<typename T, typename InputIterator, typename Predicate>
  inline
  T rms_if(InputIterator first, InputIterator last, Predicate pred)
  {
    return std::sqrt(slip::covariance_if(first,last,first,T(0),T(0),pred));
  }



  template<typename InputIterator, typename Predicate>
  inline
  typename std::iterator_traits<InputIterator>::value_type
  median_from_sorted_data_if(InputIterator first, InputIterator last, 
                                                     Predicate pred)
  {
    std::size_t mask_count=slip::cardinal_if<std::size_t>(first,last,pred);
    assert(mask_count != 0);
    std::size_t count = 0;
    while(first != last)
      {
                if (pred(*first))
                  {
                    count++;
                  }
                if (count > mask_count/2)
                  {
                    return *first;
                  }
                first++;
      }    
    std::cerr << "Error computing median_from_sorted_data_if"<<std::endl;
    return *first; 
  }
  
  template<typename InputIterator, typename Predicate>
  inline
  typename std::iterator_traits<InputIterator>::value_type
  median_if(InputIterator first, InputIterator last, 
                    Predicate pred)
  {
    std::size_t mask_count=slip::cardinal_if<std::size_t>(first,last,pred);
    assert(mask_count != 0);
    slip::Array<typename std::iterator_traits<InputIterator>::value_type> tmp(mask_count);
    std::size_t count = 0;

    while(first != last)
      {
                if (pred(*first))
                  {
                    tmp[count]=*first;
                    count++;
                  }
                first++;
      }    
    std::sort(tmp.begin(),tmp.end());
    return slip::median_from_sorted_data_n(tmp.begin(),mask_count); 

  }

  template<typename InputIterator, typename Predicate>
  inline
  typename std::iterator_traits<InputIterator>::value_type
  first_quartile_if(InputIterator first, InputIterator last, Predicate pred)
  {
    std::size_t mask_count=slip::cardinal_if<std::size_t>(first,last,pred);
    assert(mask_count != 0);
    slip::Array<typename std::iterator_traits<InputIterator>::value_type> tmp(mask_count);
    std::size_t count = 0;

    while(first != last)
      {
                if (pred(*first))
                  {
                    tmp[count] = *first;
                    count++;
                  }
                first++;
      }    

    std::sort(tmp.begin(),tmp.end());
    return slip::first_quartile(tmp.begin(),tmp.end()); 

  }


  template<typename InputIterator, typename Predicate>
  inline
  typename std::iterator_traits<InputIterator>::value_type
  third_quartile_if(InputIterator first, InputIterator last, Predicate pred)
  {
    std::size_t mask_count=slip::cardinal_if<std::size_t>(first,last,pred);
    assert(mask_count != 0);
    slip::Array<typename std::iterator_traits<InputIterator>::value_type> tmp(mask_count);
    std::size_t count = 0;

    while(first != last)
      {
                if (pred(*first))
                  {
                    tmp[count] = *first;
                    count++;
                  }
                first++;
      }    

    std::sort(tmp.begin(),tmp.end());
    return slip::third_quartile(tmp.begin(),tmp.end()); 

  }

  
  template<typename InputIterator, typename Predicate, typename OutputIterator>
  inline
  void center_if(InputIterator first, InputIterator last, 
                                 Predicate pred,
                                 OutputIterator out_first)
  {
    assert(first != last); 

    typedef typename std::iterator_traits<OutputIterator>::value_type value_type;
    value_type mean = slip::mean_if<value_type>(first,last,pred);
   
    while(first != last)
        {
                  if(pred(*first))
            { 
                      *out_first = (*first - mean);
                    }
                  ++out_first;
                  ++first;
        }
   }


  template<typename InputIterator, typename Predicate, typename OutputIterator>
  inline
  void studentize_if(InputIterator first, InputIterator last, 
                                 Predicate pred,
                                 OutputIterator out_first)
  {
    assert(first != last); 

    typedef typename std::iterator_traits<OutputIterator>::value_type value_type;
    value_type mean = slip::mean_if<value_type>(first,last,pred);
    value_type std_dev = slip::std_dev_if<value_type>(first,last,mean,pred);
   
    if(std_dev != value_type(0))
      {
                while(first != last)
                  {
                    if(pred(*first))
                      { 
                                *out_first = (*first - mean) / std_dev;
                      }
                    ++out_first;
                    ++first;
                  }
      }
    else
      {
                while(first != last)
                  {
                    if(pred(*first))
                      { 
                                *out_first = (*first - mean);
                      }
                    ++out_first;
                    ++first;
                  }
      }
   }






template <typename InputIterator, typename T>
 void statistics(InputIterator first,
                 InputIterator last,
                 slip::Statistics<T>& statistic)
 
    {   
      std::size_t n = slip::cardinal<std::size_t>(first,last);
      slip::Array<typename std::iterator_traits<InputIterator>::value_type> 
                tmp(n,first,last);
      std::sort(tmp.begin(),tmp.end());

      
      statistic.min() = *tmp.begin();
      statistic.first_quartile() = *(tmp.begin()+n/4);
      statistic.median() = *(tmp.begin()+n/2);
      statistic.third_quartile()= *((tmp.end()-1)-n/4);
      statistic.max() = *(tmp.end() - 1);
      statistic.mean() = slip::mean<T>(first,last);
      statistic.std_dev()= slip::std_dev(first,last,statistic.mean());
      statistic.skewness() = slip::skewness(first,last,statistic.mean());
      statistic.kurtosis() = slip::kurtosis(first,last,statistic.mean());
      statistic.cardinal() = n; 

     }
  

template <typename InputIterator, typename T>
 void unbiased_statistics(InputIterator first,
                                                  InputIterator last,
                                                  slip::Statistics<T>& statistic)
 
    {   
      std::size_t n = slip::cardinal<std::size_t>(first,last);
      slip::Array<typename std::iterator_traits<InputIterator>::value_type> 
                tmp(n,first,last);
      std::sort(tmp.begin(),tmp.end());

      
      statistic.min() = *tmp.begin();
      statistic.first_quartile() = *(tmp.begin()+n/4);
      statistic.median() = *(tmp.begin()+n/2);
      statistic.third_quartile()= *((tmp.end()-1)-n/4);
      statistic.max() = *(tmp.end() - 1);
      statistic.mean() = slip::mean<T>(first,last);
      statistic.std_dev()= slip::unbiased_std_dev(first,last,statistic.mean());
      statistic.skewness() = slip::unbiased_skewness(first,last,statistic.mean());
      statistic.kurtosis() = slip::unbiased_kurtosis(first,last,statistic.mean());
      statistic.cardinal() = n; 

     }
  


  template <typename InputIterator, typename T, typename MaskIterator>
 void statistics_mask(InputIterator first,
                                      InputIterator last,
                                      MaskIterator mask_first, 
                                      slip::Statistics<T>& statistic,
                                      typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
 
    {   

       std::size_t n = slip::cardinal_mask<std::size_t>(mask_first,mask_first+(last-first),value);
       assert(n != 0);
    slip::Array<typename std::iterator_traits<InputIterator>::value_type> tmp(n);
    std::size_t count = 0;
    while(first != last)
      {
                if (*mask_first == value)
                  {
                    tmp[count] = *first;
                    count++;
                  }
                first++;
        mask_first++;

      }    

    std::sort(tmp.begin(),tmp.end());
     

      
    statistic.min() = *(tmp.begin());
      statistic.first_quartile() = *(tmp.begin()+n/4);
      statistic.median() = *(tmp.begin()+n/2);
      statistic.third_quartile()= *((tmp.end()-1)-n/4);
      statistic.max() = *(tmp.end() - 1);
      statistic.mean() = slip::mean<T>(tmp.begin(),tmp.end());
      statistic.std_dev()= slip::std_dev(tmp.begin(),tmp.end(),statistic.mean());
      statistic.skewness() = slip::skewness(tmp.begin(),tmp.end(),statistic.mean());
      statistic.kurtosis() = slip::kurtosis(tmp.begin(),tmp.end(),statistic.mean());
      statistic.cardinal() = n; 

     }


  template <typename InputIterator, typename T, typename MaskIterator>
 void unbiased_statistics_mask(InputIterator first,
                                      InputIterator last,
                                      MaskIterator mask_first, 
                                      slip::Statistics<T>& statistic,
                                      typename std::iterator_traits<MaskIterator>::value_type value=typename std::iterator_traits<MaskIterator>::value_type(1))
 
    {   

       std::size_t n = slip::cardinal_mask<std::size_t>(mask_first,mask_first+(last-first),value);
       assert(n != 0);
    slip::Array<typename std::iterator_traits<InputIterator>::value_type> tmp(n);
    std::size_t count = 0;
    while(first != last)
      {
                if (*mask_first == value)
                  {
                    tmp[count] = *first;
                    count++;
                  }
                first++;
        mask_first++;

      }    

    std::sort(tmp.begin(),tmp.end());
     

      
    statistic.min() = *(tmp.begin());
      statistic.first_quartile() = *(tmp.begin()+n/4);
      statistic.median() = *(tmp.begin()+n/2);
      statistic.third_quartile()= *((tmp.end()-1)-n/4);
      statistic.max() = *(tmp.end() - 1);
      statistic.mean() = slip::mean<T>(tmp.begin(),tmp.end());
      statistic.std_dev()= slip::unbiased_std_dev(tmp.begin(),tmp.end(),statistic.mean());
      statistic.skewness() = slip::unbiased_skewness(tmp.begin(),tmp.end(),statistic.mean());
      statistic.kurtosis() = slip::unbiased_kurtosis(tmp.begin(),tmp.end(),statistic.mean());
      statistic.cardinal() = n; 

     }


  template <typename InputIterator, typename T, typename Predicate>
 void statistics_if(InputIterator first,
                                    InputIterator last,
                                    slip::Statistics<T>& statistic,
                                    Predicate pred)
                                    
    {   

       std::size_t n = slip::cardinal_if<std::size_t>(first,last,pred);
       assert(n != 0);
    slip::Array<typename std::iterator_traits<InputIterator>::value_type> tmp(n);
    std::size_t count = 0;
    while(first != last)
      {
                if (pred(*first))
                  {
                    tmp[count]=*first;
                    count++;
                }
                first++;
      }    

  
    std::sort(tmp.begin(),tmp.end());
     

      
    statistic.min() = *tmp.begin();
    statistic.first_quartile() = *(tmp.begin()+n/4);
    statistic.median() = *(tmp.begin()+n/2);
    statistic.third_quartile()= *((tmp.end()-1)-n/4);
    statistic.max() = *(tmp.end() - 1);
    statistic.mean() = slip::mean<T>(tmp.begin(),tmp.end());
    statistic.std_dev()= slip::std_dev(tmp.begin(),tmp.end(),statistic.mean());
    statistic.skewness() = slip::skewness(tmp.begin(),tmp.end(),statistic.mean());
    statistic.kurtosis() = slip::kurtosis(tmp.begin(),tmp.end(),statistic.mean());
    statistic.cardinal() = n; 

     }


  template <typename InputIterator, typename T, typename Predicate>
 void unbiased_statistics_if(InputIterator first,
                                    InputIterator last,
                                    slip::Statistics<T>& statistic,
                                    Predicate pred)
                                    
    {   

       std::size_t n = slip::cardinal_if<std::size_t>(first,last,pred);
       assert(n != 0);
    slip::Array<typename std::iterator_traits<InputIterator>::value_type> tmp(n);
    std::size_t count = 0;
    while(first != last)
      {
                if (pred(*first))
                  {
                    tmp[count]=*first;
                    count++;
                }
                first++;
      }    

  
    std::sort(tmp.begin(),tmp.end());
     

      
    statistic.min() = *tmp.begin();
    statistic.first_quartile() = *(tmp.begin()+n/4);
    statistic.median() = *(tmp.begin()+n/2);
    statistic.third_quartile()= *((tmp.end()-1)-n/4);
    statistic.max() = *(tmp.end() - 1);
    statistic.mean() = slip::mean<T>(tmp.begin(),tmp.end());
    statistic.std_dev()= slip::unbiased_std_dev(tmp.begin(),tmp.end(),statistic.mean());
    statistic.skewness() = slip::unbiased_skewness(tmp.begin(),tmp.end(),statistic.mean());
    statistic.kurtosis() = slip::unbiased_kurtosis(tmp.begin(),tmp.end(),statistic.mean());
    statistic.cardinal() = n; 

     }


  /* @{ */
  template<typename DataT,
                   typename MeanT>
  MeanT mean_next(const int n, 
                                  const MeanT& prev_mean,
                                  const DataT& xnp1)
  {
    MeanT an = static_cast<MeanT>(1)/static_cast<MeanT>(n+1);
    MeanT bn = static_cast<MeanT>(n)*an;
    return bn * prev_mean + an * xnp1;
  }
  

  template<typename DataT,
                   typename VarT>
  void var_next(const int n, 
                                const VarT& prev_mean,
                                const VarT& prev_var,
                                const DataT& xnp1,
                                VarT& next_mean,
                                VarT& next_var)
  {
    VarT an = static_cast<VarT>(1)/static_cast<VarT>(n+1);
    VarT bn = static_cast<VarT>(n)*an;
    VarT cn = bn * an;
    VarT s  = (static_cast<VarT>(xnp1) - prev_mean);
    next_mean = bn * prev_mean + an * static_cast<VarT>(xnp1);
    next_var =  bn * prev_var + cn * (s * s);
  }
/* @} */
  
}//slip::

#endif //SLIP_STATISTICS_HPP