source: trunk/yat/statistics/utility.h @ 932

Last change on this file since 932 was 932, checked in by Peter, 15 years ago

median and percentile functions now take iterators rather than vectors

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.2 KB
Line 
1#ifndef _theplu_yat_statistics_utility_
2#define _theplu_yat_statistics_utility_
3
4// $Id: utility.h 932 2007-10-05 21:03:46Z peter $
5
6/*
7  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
8  Copyright (C) 2005 Peter Johansson
9  Copyright (C) 2006 Jari Häkkinen, Markus Ringnér, Peter Johansson
10  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
11
12  This file is part of the yat library, http://trac.thep.lu.se/trac/yat
13
14  The yat library is free software; you can redistribute it and/or
15  modify it under the terms of the GNU General Public License as
16  published by the Free Software Foundation; either version 2 of the
17  License, or (at your option) any later version.
18
19  The yat library is distributed in the hope that it will be useful,
20  but WITHOUT ANY WARRANTY; without even the implied warranty of
21  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22  General Public License for more details.
23
24  You should have received a copy of the GNU General Public License
25  along with this program; if not, write to the Free Software
26  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
27  02111-1307, USA.
28*/
29
30#include "yat/classifier/DataLookupWeighted1D.h"
31#include "yat/classifier/Target.h"
32#include "yat/utility/vector.h"
33#include "yat/utility/yat_assert.h"
34
35#include <algorithm>
36#include <cmath>
37#include <vector>
38
39#include <gsl/gsl_statistics_double.h>
40
41namespace theplu {
42namespace yat {
43namespace statistics { 
44
45  //forward declarations
46  template <class T> 
47  double median(T first, T last, const bool sorted=false); 
48
49  template <class T>
50  double percentile(T first, T last, double p, bool sorted=false);
51 
52  /**
53     Adding each value in an array \a v \a to an object \a o.  The
54     requirements for the type T1 is to have an add(double, bool)
55     function, and for T2 of the array \a v are: operator[] returning
56     an element and function size() returning the number of elements.
57  */
58  template <typename T1, typename T2>
59  void add(T1& o, const T2& v, const classifier::Target& target)
60  {
61    for (size_t i=0; i<v.size(); ++i) 
62      o.add(v[i],target.binary(i));
63  } 
64
65  /**
66     Adding each value in an array \a v \a to an object \a o.  The
67     requirements for the type T1 is to have an add(double, bool)
68     function, and for T2 of the array \a v are: operator[] returning
69     an element and function size() returning the number of elements.
70  */
71  template <typename T1>
72  void add(T1& o, const classifier::DataLookupWeighted1D& v, 
73           const classifier::Target& target)
74  {
75    for (size_t i=0; i<v.size(); ++i) 
76      o.add(v.data(i),target.binary(i),v.weight(i));
77  } 
78
79  ///
80  /// Calculates the probability to get \a k or smaller from a
81  /// hypergeometric distribution with parameters \a n1 \a n2 \a
82  /// t. Hypergeomtric situation you get in the following situation:
83  /// Let there be \a n1 ways for a "good" selection and \a n2 ways
84  /// for a "bad" selection out of a total of possibilities. Take \a
85  /// t samples without replacement and \a k of those are "good"
86  /// samples. \a k will follow a hypergeomtric distribution.
87  ///
88  /// @return cumulative hypergeomtric distribution functions P(k).
89  ///
90  double cdf_hypergeometric_P(u_int k, u_int n1, u_int n2, u_int t);
91
92
93  ///
94  /// @brief Computes the kurtosis of the data in a vector.
95  ///
96  /// The kurtosis measures how sharply peaked a distribution is,
97  /// relative to its width. The kurtosis is normalized to zero for a
98  /// gaussian distribution.
99  ///
100  double kurtosis(const utility::vector&);
101
102
103  ///
104  /// @brief Median absolute deviation from median
105  ///
106  template <class T>
107  double mad(const std::vector<T>& vec, const bool sorted=false)
108  {
109    double m = median(vec, sorted);
110    std::vector<double> ad;
111    ad.reserve(vec.size());
112    for (size_t i = 0; i<vec.size(); ++i)
113      ad.push_back(fabs(vec[i]-m));
114    std::sort(ad.begin(), ad.end());
115    return median(ad.begin(), ad.end(),true);
116  }
117 
118
119  ///
120  /// @brief Median absolute deviation from median
121  ///
122  double mad(const utility::vector& vec, const bool sorted=false);
123 
124
125  ///
126  /// Median is defined to be value in the middle. If number of values
127  /// is even median is the average of the two middle values.  the
128  /// median value is given by p equal to 50. If \a sorted is false
129  /// (default), the range is copied, the copy is sorted, and then
130  /// used to calculate the median.
131  ///
132  /// Requirements: T should be an iterator over a range of doubles (or
133  /// any type being convertable to double). If \a sorted is false
134  /// iterator must be mutable, else read-only iterator is also ok.
135  ///
136  /// @return median of range
137  ///
138  template <class T> 
139  double median(T first, T last, const bool sorted=false) 
140  { return percentile(first, last, 50.0, sorted); }
141
142  /**
143     The percentile is determined by the \a p, a number between 0 and
144     100. The percentile is found by interpolation, using the formula
145     \f$ percentile = (1 - \delta) x_i + \delta x_{i+1} \f$ where \a
146     p is floor\f$((n - 1)p/100)\f$ and \f$ \delta \f$ is \f$
147     (n-1)p/100 - i \f$.Thus the minimum value of the vector is given
148     by p equal to zero, the maximum is given by p equal to 100 and
149     the median value is given by p equal to 50. If @a sorted
150     is false (default), the vector is copied, the copy is sorted,
151     and then used to calculate the median.
152
153     Requirements: T should be an iterator over a range of doubles (or
154     any type being convertable to double). If \a sorted is false
155     iterator must be mutable, else read-only iterator is also ok.
156     
157     @return \a p'th percentile of range
158  */
159  template <class T>
160  double percentile(T first, T last, double p, bool sorted=false)
161  {
162    yat_assert(first<last && "range is invalid");
163    yat_assert(p>=0);
164    yat_assert(p<=100);
165    if (sorted){
166      if (p>=100)
167        return *(--last);
168      double j = p/100 * (std::distance(first,last)-1);
169      int i = static_cast<int>(j);
170      return (1-j+floor(j))*first[i] + (j-floor(j))*first[i+1];
171    }
172
173    std::vector<double> v_copy;
174    v_copy.reserve(std::distance(first,last));
175    std::copy(first, last, std::back_inserter(v_copy));
176    double j = p/100 * (v_copy.size()-1);
177    int i = static_cast<int>(j);
178    std::partial_sort(v_copy.begin(),v_copy.begin()+i+2 , v_copy.end());
179    return percentile(v_copy.begin(), v_copy.end(), p, true);
180  }
181
182  ///
183  /// @brief Computes the skewness of the data in a vector.
184  ///
185  /// The skewness measures the asymmetry of the tails of a
186  /// distribution.
187  ///
188  double skewness(const utility::vector&);
189 
190}}} // of namespace statistics, yat, and theplu
191
192#endif
Note: See TracBrowser for help on using the repository browser.