peaks_over_threshold.hpp

来自「Boost provides free peer-reviewed portab」· HPP 代码 · 共 400 行 · 第 1/2 页

    struct peaks_over_threshold_prob_impl
      : accumulator_base
    {
        typedef typename numeric::functional::average<Sample, std::size_t>::result_type float_type;
        // for boost::result_of
        typedef boost::tuple<float_type, float_type, float_type> result_type;
        // for left tail fitting, mirror the extreme values
        typedef mpl::int_<is_same<LeftRight, left>::value ? -1 : 1> sign;

        template<typename Args>
        peaks_over_threshold_prob_impl(Args const &args)
          : mu_(sign::value * numeric::average(args[sample | Sample()], (std::size_t)1))
          , sigma2_(numeric::average(args[sample | Sample()], (std::size_t)1))
          , threshold_probability_(args[pot_threshold_probability])
          , fit_parameters_(boost::make_tuple(0., 0., 0.))
          , is_dirty_(true)
        {
        }

        void operator ()(dont_care)
        {
            this->is_dirty_ = true;
        }

        template<typename Args>
        result_type result(Args const &args) const
        {
            if (this->is_dirty_)
            {
                this->is_dirty_ = false;

                std::size_t cnt = count(args);

                // the n'th cached sample provides an approximate threshold value u
                std::size_t n = static_cast<std::size_t>(
                    std::ceil(
                        cnt * ( ( is_same<LeftRight, left>::value ) ? this->threshold_probability_ : 1. - this->threshold_probability_ )
                    )
                );

                // If n is in a valid range, return result, otherwise return NaN or throw exception
                if ( n >= static_cast<std::size_t>(tail(args).size()))
                {
                    if (std::numeric_limits<float_type>::has_quiet_NaN)
                    {
                        return boost::make_tuple(
                            std::numeric_limits<float_type>::quiet_NaN()
                          , std::numeric_limits<float_type>::quiet_NaN()
                          , std::numeric_limits<float_type>::quiet_NaN()
                        );
                    }
                    else
                    {
                        std::ostringstream msg;
                        msg << "index n = " << n << " is not in valid range [0, " << tail(args).size() << ")";
                        boost::throw_exception(std::runtime_error(msg.str()));
                        return boost::make_tuple(Sample(0), Sample(0), Sample(0));
                    }
                }
                else
                {
                    float_type u = *(tail(args).begin() + n - 1) * sign::value;

                    // compute mean and variance of samples above/under threshold value u
                    for (std::size_t i = 0; i < n; ++i)
                    {
                        mu_ += *(tail(args).begin() + i);
                        sigma2_ += *(tail(args).begin() + i) * (*(tail(args).begin() + i));
                    }

                    this->mu_ = sign::value * numeric::average(this->mu_, n);
                    this->sigma2_ = numeric::average(this->sigma2_, n);
                    this->sigma2_ -= this->mu_ * this->mu_;

                    if (is_same<LeftRight, left>::value)
                        this->threshold_probability_ = 1. - this->threshold_probability_;

                    float_type tmp = numeric::average(( this->mu_ - u )*( this->mu_ - u ), this->sigma2_);
                    float_type xi_hat = 0.5 * ( 1. - tmp );
                    float_type beta_hat = 0.5 * ( this->mu_ - u ) * ( 1. + tmp );
                    float_type beta_bar = beta_hat * std::pow(1. - threshold_probability_, xi_hat);
                    float_type u_bar = u - beta_bar * ( std::pow(1. - threshold_probability_, -xi_hat) - 1.)/xi_hat;
                    this->fit_parameters_ = boost::make_tuple(u_bar, beta_bar, xi_hat);
                }
            }

            return this->fit_parameters_;
        }

    private:
        mutable float_type mu_;                     // mean of samples above threshold u
        mutable float_type sigma2_;                 // variance of samples above threshold u
        mutable float_type threshold_probability_;
        mutable result_type fit_parameters_;        // boost::tuple that stores fit parameters
        mutable bool is_dirty_;
    };

} // namespace impl

///////////////////////////////////////////////////////////////////////////////
// tag::peaks_over_threshold
//
namespace tag
{
    template<typename LeftRight>
    struct peaks_over_threshold
      : depends_on<count>
      , pot_threshold_value
    {
        /// INTERNAL ONLY
        ///
        typedef accumulators::impl::peaks_over_threshold_impl<mpl::_1, LeftRight> impl;
    };

    template<typename LeftRight>
    struct peaks_over_threshold_prob
      : depends_on<count, tail<LeftRight> >
      , pot_threshold_probability
    {
        /// INTERNAL ONLY
        ///
        typedef accumulators::impl::peaks_over_threshold_prob_impl<mpl::_1, LeftRight> impl;
    };

    struct abstract_peaks_over_threshold
      : depends_on<>
    {
    };
}

///////////////////////////////////////////////////////////////////////////////
// extract::peaks_over_threshold
//
namespace extract
{
    extractor<tag::abstract_peaks_over_threshold> const peaks_over_threshold = {};
}

using extract::peaks_over_threshold;

// peaks_over_threshold<LeftRight>(with_threshold_value) -> peaks_over_threshold<LeftRight>
template<typename LeftRight>
struct as_feature<tag::peaks_over_threshold<LeftRight>(with_threshold_value)>
{
    typedef tag::peaks_over_threshold<LeftRight> type;
};

// peaks_over_threshold<LeftRight>(with_threshold_probability) -> peaks_over_threshold_prob<LeftRight>
template<typename LeftRight>
struct as_feature<tag::peaks_over_threshold<LeftRight>(with_threshold_probability)>
{
    typedef tag::peaks_over_threshold_prob<LeftRight> type;
};

template<typename LeftRight>
struct feature_of<tag::peaks_over_threshold<LeftRight> >
  : feature_of<tag::abstract_peaks_over_threshold>
{
};

template<typename LeftRight>
struct feature_of<tag::peaks_over_threshold_prob<LeftRight> >
  : feature_of<tag::abstract_peaks_over_threshold>
{
};

// So that peaks_over_threshold can be automatically substituted
// with weighted_peaks_over_threshold when the weight parameter is non-void.
template<typename LeftRight>
struct as_weighted_feature<tag::peaks_over_threshold<LeftRight> >
{
    typedef tag::weighted_peaks_over_threshold<LeftRight> type;
};

template<typename LeftRight>
struct feature_of<tag::weighted_peaks_over_threshold<LeftRight> >
  : feature_of<tag::peaks_over_threshold<LeftRight> >
{};

// So that peaks_over_threshold_prob can be automatically substituted
// with weighted_peaks_over_threshold_prob when the weight parameter is non-void.
template<typename LeftRight>
struct as_weighted_feature<tag::peaks_over_threshold_prob<LeftRight> >
{
    typedef tag::weighted_peaks_over_threshold_prob<LeftRight> type;
};

template<typename LeftRight>
struct feature_of<tag::weighted_peaks_over_threshold_prob<LeftRight> >
  : feature_of<tag::peaks_over_threshold_prob<LeftRight> >
{};

}} // namespace boost::accumulators

#ifdef _MSC_VER
# pragma warning(pop)
#endif

#endif