statistics.h

Go to the documentation of this file.

 
#ifndef GUL14_STATISTICS_H_
#define GUL14_STATISTICS_H_
 
#include <algorithm>
#include <cmath>
#include <limits>
#include <numeric>
#include <type_traits>
#include <vector>
 
#include "gul14/internal.h"
#include "gul14/traits.h"
 
namespace gul14 {
 
using statistics_result_type = double; 
 
template <typename ElementT>
auto ElementAccessor()
{
    return [](ElementT const& el) -> ElementT const&
        { return el; };
}
 
template <typename DataT, typename = void, typename = std::enable_if_t<std::is_arithmetic<DataT>::value>>
struct MinMax {
    DataT min{ std::numeric_limits<DataT>::max() }; 
    DataT max{ std::numeric_limits<DataT>::lowest() }; 
};
 
template <typename DataT>
struct MinMax<DataT, std::enable_if_t<std::is_floating_point<DataT>::value>> {
    DataT min{ NAN };
    DataT max{ NAN };
};
 
template <typename DataT, typename = std::enable_if_t<std::is_floating_point<DataT>::value>>
class StandardDeviationMean {
public:
    DataT sigma_{ NAN }; 
    DataT mean_{ NAN }; 
 
    operator DataT() const noexcept {
        return sigma_;
    }
    auto sigma() const noexcept -> DataT {
        return sigma_;
    }
    auto mean() const noexcept -> DataT {
        return mean_;
    }
};
 
 
template <typename ResultT = statistics_result_type,
          typename ContainerT,
          typename ElementT = typename ContainerT::value_type,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = typename std::decay_t<std::result_of_t<Accessor(ElementT)>>,
          typename = std::enable_if_t<IsContainerLike<ContainerT>::value>
         >
auto mean(ContainerT const& container, Accessor accessor = ElementAccessor<ElementT>()) -> ResultT
{
    auto const sum = std::accumulate(
            container.cbegin(), container.cend(),
            ResultT{ },
            [accessor] (ResultT const& accu, ElementT const& el) {
                return accu + static_cast<ResultT>(accessor(el)); } );
    return sum / static_cast<ResultT>(container.size());
}
 
template <typename ResultT = statistics_result_type,
          typename ContainerT,
          typename ElementT = typename ContainerT::value_type,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = typename std::decay_t<std::result_of_t<Accessor(ElementT)>>,
          typename = std::enable_if_t<IsContainerLike<ContainerT>::value>
         >
auto rms(ContainerT const& container, Accessor accessor = ElementAccessor<ElementT>()) -> ResultT
{
    auto const sum = std::accumulate(
            container.cbegin(), container.cend(),
            ResultT{ },
            [accessor] (ResultT const& accu, ElementT const& el) {
                return accu + std::pow(static_cast<ResultT>(accessor(el)), 2); } );
    return std::sqrt(sum / static_cast<ResultT>(container.size()));
}
 
template <typename ResultT = statistics_result_type,
          typename ContainerT,
          typename ElementT = typename ContainerT::value_type,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = typename std::decay_t<std::result_of_t<Accessor(ElementT)>>,
          typename = std::enable_if_t<IsContainerLike<ContainerT>::value>
         >
auto median(ContainerT const& container, Accessor accessor = ElementAccessor<ElementT>()) -> ResultT
{
    auto const len = container.size();
    if (len == 0)
        return std::numeric_limits<ResultT>::quiet_NaN();
 
    // work with a copy of the data
    // because nth_element() partially sorts the input data
    auto data_copy = std::vector<DataT>{ };
    data_copy.resize(len);
    auto data_copy_it = data_copy.begin();
    auto data_it = container.cbegin();
    auto const data_end = container.cend();
    for (; data_it != data_end; ++data_it) {
        *(data_copy_it++) = accessor(*data_it);
    }
 
    // What is the middle element?
    auto middle = data_copy.begin() + (len / 2);
    std::nth_element(data_copy.begin(), middle, data_copy.end());
    auto median = static_cast<ResultT>(*middle);
 
    // If we have an even number of elements we need to do more:
    // We calculate the mean value of the two 'middle' elements
    if (0 == len % 2) {
        std::nth_element(data_copy.begin(), middle - 1, data_copy.end());
        median = (median / 2) + (static_cast<ResultT>(*(middle - 1)) / 2);
    }
 
    return median;
}
 
template <typename ContainerT,
    typename ElementT = typename ContainerT::value_type,
    typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
    typename DataT = typename std::decay_t<std::result_of_t<Accessor(ElementT)>>,
    typename = std::enable_if_t<IsContainerLike<ContainerT>::value>
>
auto maximum(ContainerT const& container, Accessor accessor = ElementAccessor<ElementT>()) -> DataT
{
    constexpr auto initial_value = std::numeric_limits<DataT>::has_quiet_NaN ?
        std::numeric_limits<DataT>::quiet_NaN() : std::numeric_limits<DataT>::lowest();
 
    return std::accumulate(
        container.cbegin(), container.cend(), initial_value,
        [&accessor](DataT const& accu, ElementT const& el) -> DataT {
            auto const val = accessor(el);
            // Test portably for not-NAN (some compilers do not have std::isnan() for
            // integral types)
            if (val == val) {
                if (not (val <= accu)) // inverted logic to handle NAN correctly
                    return val;
            }
            return accu; });
}
 
template <typename ContainerT,
    typename ElementT = typename ContainerT::value_type,
    typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
    typename DataT = typename std::decay_t<std::result_of_t<Accessor(ElementT)>>,
    typename = std::enable_if_t<IsContainerLike<ContainerT>::value>
>
auto minimum(ContainerT const& container, Accessor accessor = ElementAccessor<ElementT>()) -> DataT
{
    constexpr auto initial_value = std::numeric_limits<DataT>::has_quiet_NaN ?
        std::numeric_limits<DataT>::quiet_NaN() : std::numeric_limits<DataT>::max();
 
    return std::accumulate(
        container.cbegin(), container.cend(), initial_value,
        [&accessor](DataT const& accu, ElementT const& el) -> DataT {
            auto const val = accessor(el);
            // Test portably for not-NAN (some compilers do not have std::isnan() for
            // integral types)
            if (val == val) {
                if (not (val >= accu)) // inverted logic to handle NAN correctly
                    return val;
            }
            return accu; });
}
 
template <typename ContainerT,
          typename ElementT = typename ContainerT::value_type,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = typename std::decay_t<std::result_of_t<Accessor(ElementT)>>,
          typename = std::enable_if_t<IsContainerLike<ContainerT>::value>
         >
auto min_max(ContainerT const& container, Accessor accessor = ElementAccessor<ElementT>()) -> MinMax<DataT>
{
    using MinMaxT = MinMax<DataT>;
    auto const sum = std::accumulate(
            container.cbegin(), container.cend(),
            MinMaxT{ },
            [accessor] (MinMaxT const& accu, ElementT const& el) -> MinMaxT {
                auto out{ accu };
                auto const val = accessor(el);
                // Test portably for not-NAN (some compilers do not have std::isnan() for
                // integral types)
                if (val == val) {
                    // (a >= NAN) and (a <= NAN) always false for all a
                    if (not (val >= out.min))
                        out.min = val;
                    if (not (val <= out.max))
                        out.max = val;
                }
                return out; });
    return sum;
}
 
template <typename ContainerT,
          typename ElementT = typename ContainerT::value_type,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = typename std::decay_t<std::result_of_t<Accessor(ElementT)>>,
          typename = std::enable_if_t<IsContainerLike<ContainerT>::value>
         >
auto remove_outliers(ContainerT&& cont, std::size_t outliers,
        Accessor accessor = ElementAccessor<ElementT>()) -> ContainerT&
{
    while (outliers-- > 0 and cont.size() > 0) {
        auto max_distant = std::max_element(cont.begin(), cont.end(),
                [mean = mean(cont, accessor), accessor] (ElementT const& a, ElementT const& b)
                { return std::abs(accessor(a) - mean) < std::abs(accessor(b) - mean); });
        cont.erase(max_distant);
    }
    return cont;
}
 
template <typename ContainerT,
          typename ElementT = typename ContainerT::value_type,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = typename std::decay_t<std::result_of_t<Accessor(ElementT)>>,
          typename = std::enable_if_t<IsContainerLike<ContainerT>::value>
         >
auto remove_outliers(ContainerT const& cont, std::size_t outliers,
        Accessor accessor = ElementAccessor<ElementT>()) -> std::vector<ElementT>
{
    auto c = std::vector<ElementT>(cont.size());
    std::copy(cont.cbegin(), cont.cend(), c.begin());
    return remove_outliers(std::move(c), outliers, accessor);
}
 
template <typename ResultT = statistics_result_type,
          typename ContainerT,
          typename ElementT = typename ContainerT::value_type,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = typename std::decay_t<std::result_of_t<Accessor(ElementT)>>,
          typename = std::enable_if_t<IsContainerLike<ContainerT>::value>
         >
auto standard_deviation(ContainerT const& container, Accessor accessor = ElementAccessor<ElementT>()) -> StandardDeviationMean<ResultT>
{
    auto const len = container.size();
 
    if (len == 0)
        return { };
 
    auto mean_val = mean<ResultT>(container, accessor);
 
    if (len == 1)
        return { std::numeric_limits<ResultT>::quiet_NaN(), mean_val };
 
    auto sum = std::accumulate(container.cbegin(), container.cend(),
        ResultT{ },
        [mean_val, accessor] (ResultT const& accu, ElementT const& el)
        { return accu + std::pow(static_cast<ResultT>(accessor(el)) - mean_val, 2); });
 
    sum /= static_cast<ResultT>(container.size() - 1);
 
    return { std::sqrt(sum), mean_val };
}
 
template <typename ResultT = statistics_result_type,
          typename ContainerT,
          typename ElementT = typename ContainerT::value_type,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = typename std::decay_t<std::result_of_t<Accessor(ElementT)>>,
          typename OpClosure,
          typename = std::enable_if_t<IsContainerLike<ContainerT>::value>
         >
auto accumulate(ContainerT const& container, OpClosure op, Accessor accessor = ElementAccessor<ElementT>()) -> ResultT
{
    auto const sum = std::accumulate(
            container.cbegin(), container.cend(),
            ResultT{ },
            [accessor, op] (ResultT const& accu, ElementT const& el) {
                return op(accu, accessor(el)); } );
    return sum;
}
 
namespace {
 
// The following stuff is only there to have a two iterator interface
 
    template <typename IteratorT>
    struct ContainerView {
        IteratorT const& begin_;
        IteratorT const& end_;
        using value_type = std::decay_t<decltype(*begin_)>;
 
        ContainerView(IteratorT const& i1, IteratorT const& i2)
        : begin_{ i1 }
        , end_{ i2 }
        {
        }
 
        // Just implement the member functions that we use here
 
        auto cbegin() const noexcept -> IteratorT const&
        {
            return begin_;
        }
        auto cend() const noexcept -> IteratorT const&
        {
            return end_;
        }
 
        auto size() const noexcept -> std::size_t
        {
            return std::distance(begin_, end_);
        }
    };
 
    template<typename IteratorT>
    auto make_view(IteratorT const& cbegin, IteratorT const& cend) -> ContainerView<IteratorT> const
    {
        return ContainerView<IteratorT>{ cbegin, cend };
    }
 
} // namespace anonymous
 
template <typename ResultT = statistics_result_type,
          typename IteratorT,
          typename ElementT = std::decay_t<decltype(*std::declval<IteratorT>())>,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = std::decay_t<std::result_of_t<Accessor(ElementT)>>>
auto mean(IteratorT const& begin, IteratorT const& end,
        Accessor accessor = ElementAccessor<ElementT>()) -> ResultT
{
    return mean<ResultT>(make_view(begin, end), accessor);
}
 
template <typename ResultT = statistics_result_type,
          typename IteratorT,
          typename ElementT = std::decay_t<decltype(*std::declval<IteratorT>())>,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = std::decay_t<std::result_of_t<Accessor(ElementT)>>>
auto rms(IteratorT const& begin, IteratorT const& end,
        Accessor accessor = ElementAccessor<ElementT>()) -> ResultT
{
    return rms<ResultT>(make_view(begin, end), accessor);
}
 
template <typename ResultT = statistics_result_type,
          typename IteratorT,
          typename ElementT = std::decay_t<decltype(*std::declval<IteratorT>())>,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = std::decay_t<std::result_of_t<Accessor(ElementT)>>>
auto median(IteratorT const& begin, IteratorT const& end,
        Accessor accessor = ElementAccessor<ElementT>()) -> ResultT
{
    return median<ResultT>(make_view(begin, end), accessor);
}
 
template <typename IteratorT,
    typename ElementT = std::decay_t<decltype(*std::declval<IteratorT>())>,
    typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
    typename DataT = std::decay_t<std::result_of_t<Accessor(ElementT)>>>
    auto maximum(IteratorT const& begin, IteratorT const& end,
        Accessor accessor = ElementAccessor<ElementT>()) -> DataT
{
    return maximum(make_view(begin, end), accessor);
}
 
template <typename IteratorT,
    typename ElementT = std::decay_t<decltype(*std::declval<IteratorT>())>,
    typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
    typename DataT = std::decay_t<std::result_of_t<Accessor(ElementT)>>>
    auto minimum(IteratorT const& begin, IteratorT const& end,
        Accessor accessor = ElementAccessor<ElementT>()) -> DataT
{
    return minimum(make_view(begin, end), accessor);
}
 
template <typename IteratorT,
          typename ElementT = std::decay_t<decltype(*std::declval<IteratorT>())>,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = std::decay_t<std::result_of_t<Accessor(ElementT)>>>
auto min_max(IteratorT const& begin, IteratorT const& end,
        Accessor accessor = ElementAccessor<ElementT>()) -> MinMax<DataT>
{
    return min_max(make_view(begin, end), accessor);
}
 
template <typename IteratorT,
          typename ElementT = std::decay_t<decltype(*std::declval<IteratorT>())>,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = std::decay_t<std::result_of_t<Accessor(ElementT)>>>
auto remove_outliers(IteratorT const& begin, IteratorT const& end,
        std::size_t outliers, Accessor accessor = ElementAccessor<ElementT>()) -> std::vector<ElementT>
{
    return remove_outliers(make_view(begin, end), outliers, accessor);
}
 
template <typename ResultT = statistics_result_type,
          typename IteratorT,
          typename ElementT = std::decay_t<decltype(*std::declval<IteratorT>())>,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = std::decay_t<std::result_of_t<Accessor(ElementT)>>>
auto standard_deviation(IteratorT const& begin, IteratorT const& end,
        Accessor accessor = ElementAccessor<ElementT>()) -> StandardDeviationMean<ResultT>
{
    return standard_deviation<ResultT>(make_view(begin, end), accessor);
}
 
template <typename ResultT = statistics_result_type,
          typename IteratorT,
          typename ElementT = std::decay_t<decltype(*std::declval<IteratorT>())>,
          typename Accessor = std::result_of_t<decltype(ElementAccessor<ElementT>())(ElementT)>(*)(ElementT const&),
          typename DataT = std::decay_t<std::result_of_t<Accessor(ElementT)>>,
          typename OpClosure>
auto accumulate(IteratorT const& begin, IteratorT const& end, OpClosure op,
        Accessor accessor = ElementAccessor<ElementT>()) -> ResultT
{
    return accumulate<ResultT>(make_view(begin, end), op, accessor);
}
 
} // namespace gul14
 
#endif
 
// vi:ts=4:sw=4:sts=4:et