TukeyBiweightEstimator.h

#ifndef _theplu_yat_statistics_tukey_biweight_estimator 
#define _theplu_yat_statistics_tukey_biweight_estimator

// $Id: TukeyBiweightEstimator.h 2511 2011-07-08 23:44:33Z peter $

/*
  Copyright (C) 2011 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include "AveragerWeighted.h"
#include "utility.h"

#include "yat/regression/TukeyBiweight.h"

#include "yat/utility/concept_check.h"
#include "yat/utility/Exception.h"
#include "yat/utility/iterator_traits.h"

#include <boost/concept_check.hpp>

#include <algorithm>
#include <iterator>
#include <limits>
#include <vector>

namespace theplu {
namespace yat {
namespace statistics {  

  class TukeyBiweightEstimator
  {
  public:
    explicit TukeyBiweightEstimator(double cutoff=4.685, bool sorted=false)
      : cutoff_(cutoff), sorted_(sorted) {}

    template<typename RandomAccessIterator>
    double operator()(RandomAccessIterator first, 
                      RandomAccessIterator last) const;

  private:
    double cutoff_;
    bool sorted_;

    template<typename RandomAccessIterator>
    double estimate(RandomAccessIterator first, 
                    RandomAccessIterator last) const;

    template<typename InputIterator>
    double estimate(InputIterator first, InputIterator last, 
                    double center, double spread) const;
  };

  template<typename RandomAccessIterator>
  double TukeyBiweightEstimator::operator()(RandomAccessIterator first, 
                                            RandomAccessIterator last) const
  {
    BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<RandomAccessIterator>));
    using utility::DataIteratorConcept;
    BOOST_CONCEPT_ASSERT((DataIteratorConcept<RandomAccessIterator>));

    if (sorted_)
      return estimate(first, last);

    // if not sorted, create a sorted copy
    typedef typename std::iterator_traits<RandomAccessIterator> traits;
    std::vector<typename traits::value_type> vec;
    vec.reserve(last-first);
    std::copy(first, last, std::back_inserter(vec));
    std::sort(vec.begin(), vec.end());
    return estimate(vec.begin(), vec.end());

    const double scale = mad(first, last, true);
    double m = median(first, last, true);
    if (scale==0)
      return m;
    // ensure m != m0
    double m0 = m+1.0;
    size_t epoch = 0;
    while (m!=m0) {
      m0 = m;
      m = estimate(first, last, m, scale);
      ++epoch;
      if (epoch>1000)
        utility::runtime_error("TukeyBiweightIterator: too many epochs");
    }
    return m;
  }


  template<typename RandomAccessIterator>
  double TukeyBiweightEstimator::estimate(RandomAccessIterator first, 
                                          RandomAccessIterator last) const
  {
    BOOST_CONCEPT_ASSERT((boost::RandomAccessIterator<RandomAccessIterator>));
    using utility::DataIteratorConcept;
    BOOST_CONCEPT_ASSERT((DataIteratorConcept<RandomAccessIterator>));

    const double scale = mad(first, last, true);
    double m = median(first, last, true);
    // if mad is zero all (non-zero weight) data points are equal and
    // median is the only sensible estimate. Also, calculations below
    // would "divide by zero" so we need to interupt here.
    if (scale==0)
      return m;
    double m0 = m+1.0;
    size_t epoch = 0;
    // FIXME: let user define convergence requirement
    while (m!=m0) {
      m0 = m;
      m = estimate(first, last, m, scale);
      ++epoch;
      // FIXME: let user define maximal epochs
      if (epoch>1000)
        utility::runtime_error("TukeyBiweightIterator: too many epochs");
    }
    return m;
  }


  template<typename InputIterator>
  double TukeyBiweightEstimator::estimate(InputIterator first, 
                                          InputIterator last, 
                                          double center, double spread) const
  {
    double scale = spread*cutoff_;
    AveragerWeighted averager;
    regression::TukeyBiweight biweight;
    utility::iterator_traits<InputIterator> traits;
    for ( ; first!=last; ++first) {
      double x = traits.data(first);
      double w = traits.weight(first) * biweight((x-center)/scale);
      averager.add(x, w); 
    }
    return averager.mean();
  }

}}} // end of namespace theplu yat statistics
# endif