yat  0.14.5pre
TukeyBiweightEstimator.h
1 #ifndef _theplu_yat_statistics_tukey_biweight_estimator
2 #define _theplu_yat_statistics_tukey_biweight_estimator
3 
4 // $Id: TukeyBiweightEstimator.h 3550 2017-01-03 05:41:02Z peter $
5 
6 /*
7  Copyright (C) 2011, 2012, 2014, 2016 Peter Johansson
8 
9  This file is part of the yat library, http://dev.thep.lu.se/yat
10 
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 3 of the
14  License, or (at your option) any later version.
15 
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with yat. If not, see <http://www.gnu.org/licenses/>.
23 */
24 
25 #include "AveragerWeighted.h"
26 #include "utility.h"
27 
28 #include "yat/regression/TukeyBiweight.h"
29 
30 #include "yat/utility/concept_check.h"
31 #include "yat/utility/Exception.h"
32 #include "yat/utility/iterator_traits.h"
33 
34 #include <boost/concept_check.hpp>
35 #include <boost/iterator/iterator_concepts.hpp>
36 
37 #include <algorithm>
38 #include <iterator>
39 #include <limits>
40 #include <vector>
41 
42 namespace theplu {
43 namespace yat {
44 namespace statistics {
45 
65  {
66  public:
76  explicit TukeyBiweightEstimator(double cutoff=4.685, bool sorted=false)
77  : cutoff_(cutoff), sorted_(sorted) {}
78 
87  template<typename RandomAccessIterator>
88  double operator()(RandomAccessIterator first,
89  RandomAccessIterator last) const;
90 
91  private:
92  double cutoff_;
93  bool sorted_;
94 
95  template<typename RandomAccessIterator>
96  double estimate(RandomAccessIterator first,
97  RandomAccessIterator last) const;
98 
99  template<typename InputIterator>
100  double estimate(InputIterator first, InputIterator last,
101  double center, double spread) const;
102  };
103 
104  template<typename RandomAccessIterator>
105  double TukeyBiweightEstimator::operator()(RandomAccessIterator first,
106  RandomAccessIterator last) const
107  {
108  using boost_concepts::RandomAccessTraversal;
109  BOOST_CONCEPT_ASSERT((RandomAccessTraversal<RandomAccessIterator>));
111  BOOST_CONCEPT_ASSERT((DataIteratorConcept<RandomAccessIterator>));
112 
113  if (sorted_)
114  return estimate(first, last);
115 
116  // if not sorted, create a sorted copy
117  typedef typename std::iterator_traits<RandomAccessIterator> traits;
118  std::vector<typename traits::value_type> vec(first, last);
119  std::sort(vec.begin(), vec.end());
120  return estimate(vec.begin(), vec.end());
121  }
122 
123 
124  template<typename RandomAccessIterator>
125  double TukeyBiweightEstimator::estimate(RandomAccessIterator first,
126  RandomAccessIterator last) const
127  {
128  const double scale = mad(first, last, true);
129  double m = median(first, last, true);
130  // if mad is zero all (non-zero weight) data points are equal and
131  // median is the only sensible estimate. Also, calculations below
132  // would "divide by zero" so we need to interupt here.
133  if (scale==0)
134  return m;
135  double m0 = m+1.0;
136  size_t epoch = 0;
137  // FIXME: let user define convergence requirement
138  while (m!=m0) {
139  m0 = m;
140  m = estimate(first, last, m, scale);
141  ++epoch;
142  // FIXME: let user define maximal epochs
143  if (epoch>1000)
144  utility::runtime_error("TukeyBiweightIterator: too many epochs");
145  }
146  return m;
147  }
148 
149 
150  template<typename InputIterator>
151  double TukeyBiweightEstimator::estimate(InputIterator first,
152  InputIterator last,
153  double center, double spread) const
154  {
155  double scale = spread*cutoff_;
156  AveragerWeighted averager;
157  regression::TukeyBiweight biweight;
158  utility::iterator_traits<InputIterator> traits;
159  for ( ; first!=last; ++first) {
160  double x = traits.data(first);
161  double w = traits.weight(first) * biweight((x-center)/scale);
162  averager.add(x, w);
163  }
164  return averager.mean();
165  }
166 
167 }}} // end of namespace theplu yat statistics
168 # endif
TukeyBiweightEstimator(double cutoff=4.685, bool sorted=false)
Constructor.
Definition: TukeyBiweightEstimator.h:76
double mad(RandomAccessIterator first, RandomAccessIterator last, bool sorted=false)
Median absolute deviation from median.
Definition: utility.h:390
double operator()(RandomAccessIterator first, RandomAccessIterator last) const
Definition: TukeyBiweightEstimator.h:105
Concept check for Data Iterator.
Definition: concept_check.h:228
double median(RandomAccessIterator first, RandomAccessIterator last, bool sorted=false)
Definition: utility.h:413
Class used for all runtime error detected within yat library.
Definition: Exception.h:38
Tukey&#39;s Biweight Estimator.
Definition: TukeyBiweightEstimator.h:64

Generated on Tue Sep 26 2017 02:33:29 for yat by  doxygen 1.8.5