source: trunk/lib/statistics/utility.h @ 519

// $Id: utility.h 519 2006-02-22 09:28:12Z peter $

#ifndef _theplu_statistics_utility_ 
#define _theplu_statistics_utility_ 

#include <c++_tools/gslapi/vector.h>

#include <algorithm>
#include <cassert>
#include <cmath>
#include <vector>

namespace theplu {
namespace statistics {  

  //forward declarations
  template <class T>
  double percentile(const std::vector<T>& vec, const double p, 
                    const bool sorted=false);
  

  ///
  /// Calculates the probabilty to get \a k or smaller from a
  /// hypergeometric distribution with parameters \a n1 \a n2 \a
  /// t. Hypergeomtric situation you get in the following situation:
  /// Let there be \a n1 ways for a "good" selection and \a n2 ways
  /// for a "bad" selection out of a total of possibilities. Take \a
  /// t samples without replacement and \a k of those are "good"
  /// samples. \a k will follow a hypergeomtric distribution.
  /// @cumulative hypergeomtric distribution functions P(k).
  ///
  double cdf_hypergeometric_P(u_int k, u_int n1, u_int n2, u_int t);


  ///
  /// Median is defined to be value in the middle. If number of values
  /// is even median is the average of the two middle values.  the
  /// median value is given by p equal to 50. If @a sorted is false
  /// (default), the vector is copied, the copy is sorted, and then
  /// used to calculate the median.
  ///
  /// @return median
  ///
  /// @note interface will change
  ///
  template <class T> 
  inline double median(const std::vector<T>& v, const bool sorted=false) 
  { return percentile(v, 50.0, sorted); }

  ///
  /// Median is defined to be value in the middle. If number of values
  /// is even median is the average of the two middle values. If @a
  /// sorted is true, the function assumes vector @a vec to be
  /// sorted. If @a sorted is false, the vector is copied, the copy is
  /// sorted (default), and then used to calculate the median.
  ///
  /// @return median
  ///
  double median(const gslapi::vector& vec, const bool sorted=false);

  ///
  /// The percentile is determined by the \a p, a number between 0 and
  /// 100. The percentile is found by interpolation, using the formula
  /// \f$ percentile = (1 - \delta) x_i + \delta x_{i+1} \f$ where \a
  /// p is floor\f$((n - 1)p/100)\f$ and \f$ \delta \f$ is \f$
  /// (n-1)p/100 - i \f$.Thus the minimum value of the vector is given
  /// by p equal to zero, the maximum is given by p equal to 100 and
  /// the median value is given by p equal to 50. If @a sorted
  /// is false (default), the vector is copied, the copy is sorted,
  /// and then used to calculate the median.
  /// 
  /// @return \a p'th percentile
  ///
  template <class T>
  double percentile(const std::vector<T>& vec, const double p, 
                    const bool sorted=false)
  {
    assert(!(p>100 && p<0));
    if (sorted){
      if (p>=100)
        return vec.back();
      double j = p/100 * (vec.size()-1);
      int i = static_cast<int>(j);
      return (1-j+floor(j))*vec[i] + (j-floor(j))*vec[i+1];
    }
    if (p==100)
      return  *std::max_element(vec.begin(),vec.end());
    std::vector<T> v_copy(vec);
    double j = p/100 * (v_copy.size()-1);
    int i = static_cast<int>(j);
    std::partial_sort(v_copy.begin(),v_copy.begin()+i+2 , v_copy.end());
    return (1-j+floor(j))*v_copy[i] + (j-floor(j))*v_copy[i+1];
  
  }

  ///
  /// The percentile is determined by the \a p, a number between 0 and
  /// 100. The percentile is found by interpolation, using the formula
  /// \f$ percentile = (1 - \delta) x_i + \delta x_{i+1} \f$ where \a
  /// p is floor\f$((n - 1)p/100)\f$ and \f$ \delta \f$ is \f$
  /// (n-1)p/100 - i \f$.Thus the minimum value of the vector is given
  /// by p equal to zero, the maximum is given by p equal to 100 and
  /// the median value is given by p equal to 50. If @a sorted
  /// is false (default), the vector is copied, the copy is sorted,
  /// and then used to calculate the median.
  /// 
  /// @return \a p'th percentile
  ///
  double percentile(const gslapi::vector& vec, const double, 
                    const bool sorted=false);

}} // of namespace statistics and namespace theplu

#endif