TukeyBiweightEstimator.h

#ifndef _theplu_yat_statistics_tukey_biweight_estimator
#define _theplu_yat_statistics_tukey_biweight_estimator

// $Id: TukeyBiweightEstimator.h 3925 2020-07-05 23:00:55Z peter $

/*
  Copyright (C) 2011, 2012, 2014, 2016, 2020 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 <boost/iterator/iterator_concepts.hpp>

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

namespace theplu {
namespace yat {
namespace statistics {

  class TukeyBiweightEstimator
  {
  public:
    class Stopper
    {
    public:
      Stopper(size_t max_epochs=1000);
      virtual bool operator()(double current, double previous) const=0;
    protected:
      size_t epochs(void) const;
      size_t max_epochs(void) const;
    private:
      size_t epochs_;
      size_t max_epochs_;
      friend class TukeyBiweightEstimator;
      bool stop(double current, double previous);
    };

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

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

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

  private:
    double cutoff_;
    bool sorted_;

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

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

    class DefaultStopper : public Stopper
    {
    public:
      bool operator()(double current, double previous) const;
    };
  };


  template<typename RandomAccessIterator>
  double TukeyBiweightEstimator::operator()(RandomAccessIterator first,
                                            RandomAccessIterator last) const
  {
    TukeyBiweightEstimator::DefaultStopper stopper;
    return (*this)(first, last, stopper);
  }


  template<typename RandomAccessIterator>
  double
  TukeyBiweightEstimator::operator()(RandomAccessIterator first,
                                     RandomAccessIterator last,
                                     TukeyBiweightEstimator::Stopper& stopper) const
  {
    using boost_concepts::RandomAccessTraversal;
    BOOST_CONCEPT_ASSERT((RandomAccessTraversal<RandomAccessIterator>));
    using utility::DataIteratorConcept;
    BOOST_CONCEPT_ASSERT((DataIteratorConcept<RandomAccessIterator>));

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

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


  template<typename RandomAccessIterator>
  double
  TukeyBiweightEstimator::estimate(RandomAccessIterator first,
                                   RandomAccessIterator last,
                                   TukeyBiweightEstimator::Stopper& stopper) const
  {
    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;
    while (!stopper.stop(m, m0)) {
      m0 = m;
      m = estimate(first, last, m, scale);
    }
    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