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

Last change on this file since 1025 was 1025, checked in by Peter, 14 years ago

passing VectorBase? in regression::utility - refs #256

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.6 KB
Line 
1#ifndef _theplu_yat_statistics_utility_
2#define _theplu_yat_statistics_utility_
3
4// $Id: utility.h 1025 2008-02-01 18:25:35Z 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/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/VectorBase.h"
33#include "yat/utility/yat_assert.h"
34
35#include <algorithm>
36#include <cmath>
37#include <stdexcept>
38#include <vector>
39
40#include <gsl/gsl_statistics_double.h>
41
42namespace theplu {
43namespace yat {
44namespace statistics { 
45
46  //forward declarations
47  template <class T> 
48  double median(T first, T last, const bool sorted=false); 
49
50  template <class T>
51  double percentile(T first, T last, double p, bool sorted=false);
52 
53  /**
54     Adding each value in an array \a v \a to an object \a o.  The
55     requirements for the type T1 is to have an add(double, bool)
56     function, and for T2 of the array \a v are: operator[] returning
57     an element and function size() returning the number of elements.
58  */
59  template <typename T1, typename T2>
60  void add(T1& o, const T2& v, const classifier::Target& target)
61  {
62    for (size_t i=0; i<v.size(); ++i) 
63      o.add(v[i],target.binary(i));
64  } 
65
66  /**
67     Adding each value in an array \a v \a to an object \a o.  The
68     requirements for the type T1 is to have an add(double, bool)
69     function, and for T2 of the array \a v are: operator[] returning
70     an element and function size() returning the number of elements.
71  */
72  template <typename T1>
73  void add(T1& o, const classifier::DataLookupWeighted1D& v, 
74           const classifier::Target& target)
75  {
76    for (size_t i=0; i<v.size(); ++i) 
77      o.add(v.data(i),target.binary(i),v.weight(i));
78  } 
79
80  ///
81  /// Calculates the probability to get \a k or smaller from a
82  /// hypergeometric distribution with parameters \a n1 \a n2 \a
83  /// t. Hypergeomtric situation you get in the following situation:
84  /// Let there be \a n1 ways for a "good" selection and \a n2 ways
85  /// for a "bad" selection out of a total of possibilities. Take \a
86  /// t samples without replacement and \a k of those are "good"
87  /// samples. \a k will follow a hypergeomtric distribution.
88  ///
89  /// @return cumulative hypergeomtric distribution functions P(k).
90  ///
91  double cdf_hypergeometric_P(u_int k, u_int n1, u_int n2, u_int t);
92
93
94  ///
95  /// @brief Computes the kurtosis of the data in a vector.
96  ///
97  /// The kurtosis measures how sharply peaked a distribution is,
98  /// relative to its width. The kurtosis is normalized to zero for a
99  /// gaussian distribution.
100  ///
101  double kurtosis(const utility::VectorBase&);
102
103
104  ///
105  /// @brief Median absolute deviation from median
106  ///
107  /// Function is non-mutable function
108  ///
109  template <class T>
110  double mad(T first, T last, const bool sorted=false)
111  {
112    double m = median(first, last, sorted);
113    std::vector<double> ad;
114    ad.reserve(std::distance(first, last));
115    for( ; first!=last; ++first)
116      ad.push_back(fabs(*first-m));
117    std::sort(ad.begin(), ad.end());
118    return median(ad.begin(), ad.end(), true);
119  }
120 
121
122  ///
123  /// Median is defined to be value in the middle. If number of values
124  /// is even median is the average of the two middle values.  the
125  /// median value is given by p equal to 50. If \a sorted is false
126  /// (default), the range is copied, the copy is sorted, and then
127  /// used to calculate the median.
128  ///
129  /// Function is a non-mutable function, i.e., \a first and \a last
130  /// can be const_iterators.
131  ///
132  /// Requirements: T should be an iterator over a range of doubles (or
133  /// any type being convertable to double).
134  ///
135  /// @return median of range
136  ///
137  template <class T> 
138  double median(T first, T last, const bool sorted=false) 
139  { return percentile(first, last, 50.0, sorted); }
140
141  /**
142     The percentile is determined by the \a p, a number between 0 and
143     100. The percentile is found by interpolation, using the formula
144     \f$ percentile = (1 - \delta) x_i + \delta x_{i+1} \f$ where \a
145     p is floor\f$((n - 1)p/100)\f$ and \f$ \delta \f$ is \f$
146     (n-1)p/100 - i \f$.Thus the minimum value of the vector is given
147     by p equal to zero, the maximum is given by p equal to 100 and
148     the median value is given by p equal to 50. If @a sorted
149     is false (default), the vector is copied, the copy is sorted,
150     and then used to calculate the median.
151
152     Function is a non-mutable function, i.e., \a first and \a last
153     can be const_iterators.
154     
155     Requirements: T should be an iterator over a range of doubles (or
156     any type being convertable to double). If \a sorted is false
157     iterator must be mutable, else read-only iterator is also ok.
158     
159     @return \a p'th percentile of range
160  */
161  template <class T>
162  double percentile(T first, T last, double p, bool sorted=false)
163  {
164    utility::yat_assert<std::range_error>(first<last);
165    utility::yat_assert<std::runtime_error>(p>=0, "percentage is negative");
166    utility::yat_assert<std::runtime_error>(p<=100, 
167                                            "percentage is larger than 100");
168    if (sorted){
169      if (p>=100)
170        return *(--last);
171      // range is one value only is a special case
172      if (std::distance(first, last)==1) 
173        return *first;
174      double j = p/100 * (std::distance(first,last)-1);
175      int i = static_cast<int>(j);
176      return (1-j+floor(j))*first[i] + (j-floor(j))*first[i+1];
177    }
178
179    std::vector<double> v_copy;
180    v_copy.reserve(std::distance(first,last));
181    std::copy(first, last, std::back_inserter(v_copy));
182    double j = p/100 * (v_copy.size()-1);
183    int i = static_cast<int>(j);
184    std::partial_sort(v_copy.begin(),v_copy.begin()+i+2 , v_copy.end());
185    return percentile(v_copy.begin(), v_copy.end(), p, true);
186  }
187
188  ///
189  /// @brief Computes the skewness of the data in a vector.
190  ///
191  /// The skewness measures the asymmetry of the tails of a
192  /// distribution.
193  ///
194  double skewness(const utility::VectorBase&);
195 
196}}} // of namespace statistics, yat, and theplu
197
198#endif
Note: See TracBrowser for help on using the repository browser.