basic_statistics.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <functional>
#include <iterator>
#include <limits>
#include <map>
#include <numeric>
#include <stdexcept>
#include <type_traits>
#include <vector>
 
namespace statcpp {
 
// ============================================================================
// Sum
// ============================================================================
 
template <typename Iterator>
auto sum(Iterator first, Iterator last)
{
    using value_type = typename std::iterator_traits<Iterator>::value_type;
    return std::accumulate(first, last, value_type{});
}
 
template <typename Iterator, typename Projection>
auto sum(Iterator first, Iterator last, Projection proj)
{
    using result_type = std::invoke_result_t<Projection,
        typename std::iterator_traits<Iterator>::value_type>;
    result_type total{};
    for (auto it = first; it != last; ++it) {
        total += std::invoke(proj, *it);
    }
    return total;
}
 
// ============================================================================
// Count
// ============================================================================
 
template <typename Iterator>
std::size_t count(Iterator first, Iterator last)
{
    return static_cast<std::size_t>(std::distance(first, last));
}
 
// ============================================================================
// Mean (Arithmetic Mean)
// ============================================================================
 
template <typename Iterator>
double mean(Iterator first, Iterator last)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::mean: empty range");
    }
    double total = 0.0;
    for (auto it = first; it != last; ++it) {
        total += static_cast<double>(*it);
    }
    return total / static_cast<double>(n);
}
 
template <typename Iterator, typename Projection>
double mean(Iterator first, Iterator last, Projection proj)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::mean: empty range");
    }
    double total = 0.0;
    for (auto it = first; it != last; ++it) {
        total += static_cast<double>(std::invoke(proj, *it));
    }
    return total / static_cast<double>(n);
}
 
// ============================================================================
// Median
// ============================================================================
 
template <typename Iterator>
double median(Iterator first, Iterator last)
{
    static_assert(
        std::is_base_of_v<
            std::random_access_iterator_tag,
            typename std::iterator_traits<Iterator>::iterator_category>,
        "statcpp::median requires random access iterators");
 
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::median: empty range");
    }
 
    auto mid = n / 2;
    if (n % 2 == 0) {
        return (static_cast<double>(*(first + (mid - 1))) + static_cast<double>(*(first + mid))) / 2.0;
    }
    return static_cast<double>(*(first + mid));
}
 
template <typename Iterator, typename Projection>
double median(Iterator first, Iterator last, Projection proj)
{
    static_assert(
        std::is_base_of_v<
            std::random_access_iterator_tag,
            typename std::iterator_traits<Iterator>::iterator_category>,
        "statcpp::median requires random access iterators");
 
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::median: empty range");
    }
 
    auto mid = n / 2;
    if (n % 2 == 0) {
        return (static_cast<double>(std::invoke(proj, *(first + (mid - 1))))
              + static_cast<double>(std::invoke(proj, *(first + mid)))) / 2.0;
    }
    return static_cast<double>(std::invoke(proj, *(first + mid)));
}
 
// ============================================================================
// Mode
// ============================================================================
 
template <typename Iterator>
auto mode(Iterator first, Iterator last)
{
    using value_type = typename std::iterator_traits<Iterator>::value_type;
    if (first == last) {
        throw std::invalid_argument("statcpp::mode: empty range");
    }
 
    std::map<value_type, std::size_t> freq;
    for (auto it = first; it != last; ++it) {
        // Skip NaN values (floating-point types only)
        if constexpr (std::is_floating_point_v<value_type>) {
            if (std::isnan(*it)) continue;
        }
        ++freq[*it];
    }
 
    if (freq.empty()) {
        throw std::invalid_argument("statcpp::mode: all values are NaN");
    }
 
    auto best = freq.begin();
    for (auto it = freq.begin(); it != freq.end(); ++it) {
        if (it->second > best->second) {
            best = it;
        }
    }
    return best->first;
}
 
template <typename Iterator, typename Projection>
auto mode(Iterator first, Iterator last, Projection proj)
{
    using result_type = std::invoke_result_t<Projection,
        typename std::iterator_traits<Iterator>::value_type>;
    if (first == last) {
        throw std::invalid_argument("statcpp::mode: empty range");
    }
 
    std::map<result_type, std::size_t> freq;
    for (auto it = first; it != last; ++it) {
        // Skip NaN values (floating-point types only)
        auto val = std::invoke(proj, *it);
        if constexpr (std::is_floating_point_v<result_type>) {
            if (std::isnan(val)) continue;
        }
        ++freq[val];
    }
 
    if (freq.empty()) {
        throw std::invalid_argument("statcpp::mode: all values are NaN");
    }
 
    auto best = freq.begin();
    for (auto it = freq.begin(); it != freq.end(); ++it) {
        if (it->second > best->second) {
            best = it;
        }
    }
    return best->first;
}
 
// ============================================================================
// Modes (multiple modes)
// ============================================================================
 
template <typename Iterator>
auto modes(Iterator first, Iterator last)
    -> std::vector<typename std::iterator_traits<Iterator>::value_type>
{
    using value_type = typename std::iterator_traits<Iterator>::value_type;
    if (first == last) {
        throw std::invalid_argument("statcpp::modes: empty range");
    }
 
    std::map<value_type, std::size_t> freq;
    for (auto it = first; it != last; ++it) {
        // Skip NaN values (floating-point types only)
        if constexpr (std::is_floating_point_v<value_type>) {
            if (std::isnan(*it)) continue;
        }
        ++freq[*it];
    }
 
    if (freq.empty()) {
        throw std::invalid_argument("statcpp::modes: all values are NaN");
    }
 
    // Find maximum frequency
    std::size_t max_freq = 0;
    for (const auto& pair : freq) {
        if (pair.second > max_freq) {
            max_freq = pair.second;
        }
    }
 
    // Collect all values with maximum frequency (ascending order due to map)
    std::vector<value_type> result;
    for (const auto& pair : freq) {
        if (pair.second == max_freq) {
            result.push_back(pair.first);
        }
    }
 
    return result;
}
 
template <typename Iterator, typename Projection>
auto modes(Iterator first, Iterator last, Projection proj)
    -> std::vector<std::invoke_result_t<Projection,
        typename std::iterator_traits<Iterator>::value_type>>
{
    using result_type = std::invoke_result_t<Projection,
        typename std::iterator_traits<Iterator>::value_type>;
    if (first == last) {
        throw std::invalid_argument("statcpp::modes: empty range");
    }
 
    std::map<result_type, std::size_t> freq;
    for (auto it = first; it != last; ++it) {
        // Skip NaN values (floating-point types only)
        auto val = std::invoke(proj, *it);
        if constexpr (std::is_floating_point_v<result_type>) {
            if (std::isnan(val)) continue;
        }
        ++freq[val];
    }
 
    if (freq.empty()) {
        throw std::invalid_argument("statcpp::modes: all values are NaN");
    }
 
    // Find maximum frequency
    std::size_t max_freq = 0;
    for (const auto& pair : freq) {
        if (pair.second > max_freq) {
            max_freq = pair.second;
        }
    }
 
    // Collect all values with maximum frequency (ascending order due to map)
    std::vector<result_type> result;
    for (const auto& pair : freq) {
        if (pair.second == max_freq) {
            result.push_back(pair.first);
        }
    }
 
    return result;
}
 
// ============================================================================
// Geometric Mean
// ============================================================================
 
template <typename Iterator>
double geometric_mean(Iterator first, Iterator last)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::geometric_mean: empty range");
    }
    double log_sum = 0.0;
    for (auto it = first; it != last; ++it) {
        if (static_cast<double>(*it) <= 0.0) {
            throw std::invalid_argument("statcpp::geometric_mean: all values must be positive");
        }
        log_sum += std::log(static_cast<double>(*it));
    }
    return std::exp(log_sum / static_cast<double>(n));
}
 
template <typename Iterator, typename Projection>
double geometric_mean(Iterator first, Iterator last, Projection proj)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::geometric_mean: empty range");
    }
    double log_sum = 0.0;
    for (auto it = first; it != last; ++it) {
        double val = static_cast<double>(std::invoke(proj, *it));
        if (val <= 0.0) {
            throw std::invalid_argument("statcpp::geometric_mean: all values must be positive");
        }
        log_sum += std::log(val);
    }
    return std::exp(log_sum / static_cast<double>(n));
}
 
// ============================================================================
// Harmonic Mean
// ============================================================================
 
template <typename Iterator>
double harmonic_mean(Iterator first, Iterator last)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::harmonic_mean: empty range");
    }
    double reciprocal_sum = 0.0;
    for (auto it = first; it != last; ++it) {
        double val = static_cast<double>(*it);
        // Reject both exact zero and values so small that 1/val overflows to infinity.
        // std::numeric_limits<double>::min() is the smallest positive normal double (~2.2e-308),
        // so checking |val| < min() catches subnormals and zero without false positives on
        // legitimate small values that are still representable.
        if (std::abs(val) < std::numeric_limits<double>::min()) {
            throw std::invalid_argument("statcpp::harmonic_mean: zero or near-zero value encountered");
        }
        reciprocal_sum += 1.0 / val;
    }
    return static_cast<double>(n) / reciprocal_sum;
}
 
template <typename Iterator, typename Projection>
double harmonic_mean(Iterator first, Iterator last, Projection proj)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::harmonic_mean: empty range");
    }
    double reciprocal_sum = 0.0;
    for (auto it = first; it != last; ++it) {
        double val = static_cast<double>(std::invoke(proj, *it));
        if (std::abs(val) < std::numeric_limits<double>::min()) {
            throw std::invalid_argument("statcpp::harmonic_mean: zero or near-zero value encountered");
        }
        reciprocal_sum += 1.0 / val;
    }
    return static_cast<double>(n) / reciprocal_sum;
}
 
// ============================================================================
// Trimmed Mean
// ============================================================================
 
template <typename Iterator>
double trimmed_mean(Iterator first, Iterator last, double proportion)
{
    static_assert(
        std::is_base_of_v<
            std::random_access_iterator_tag,
            typename std::iterator_traits<Iterator>::iterator_category>,
        "statcpp::trimmed_mean requires random access iterators");
 
    if (proportion < 0.0 || proportion >= 0.5) {
        throw std::invalid_argument("statcpp::trimmed_mean: proportion must be in [0.0, 0.5)");
    }
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::trimmed_mean: empty range");
    }
 
    auto trim_count = static_cast<std::size_t>(static_cast<double>(n) * proportion);
    if (n - 2 * trim_count == 0) {
        throw std::invalid_argument("statcpp::trimmed_mean: all elements trimmed");
    }
 
    double total = 0.0;
    for (std::size_t i = trim_count; i < n - trim_count; ++i) {
        total += static_cast<double>(*(first + i));
    }
    return total / static_cast<double>(n - 2 * trim_count);
}
 
template <typename Iterator, typename Projection>
double trimmed_mean(Iterator first, Iterator last, double proportion, Projection proj)
{
    static_assert(
        std::is_base_of_v<
            std::random_access_iterator_tag,
            typename std::iterator_traits<Iterator>::iterator_category>,
        "statcpp::trimmed_mean requires random access iterators");
 
    if (proportion < 0.0 || proportion >= 0.5) {
        throw std::invalid_argument("statcpp::trimmed_mean: proportion must be in [0.0, 0.5)");
    }
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::trimmed_mean: empty range");
    }
 
    auto trim_count = static_cast<std::size_t>(static_cast<double>(n) * proportion);
    if (n - 2 * trim_count == 0) {
        throw std::invalid_argument("statcpp::trimmed_mean: all elements trimmed");
    }
 
    double total = 0.0;
    for (std::size_t i = trim_count; i < n - trim_count; ++i) {
        total += static_cast<double>(std::invoke(proj, *(first + i)));
    }
    return total / static_cast<double>(n - 2 * trim_count);
}
 
// ============================================================================
// Weighted Mean
// ============================================================================
 
template <typename Iterator, typename WeightIterator>
double weighted_mean(Iterator first, Iterator last, WeightIterator weight_first, WeightIterator weight_last)
{
    if (std::distance(first, last) != std::distance(weight_first, weight_last)) {
        throw std::invalid_argument("statcpp::weighted_mean: data and weight ranges must have the same size");
    }
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_mean: empty range");
    }
 
    double sum_weighted = 0.0;
    double sum_weights = 0.0;
    auto weight_it = weight_first;
 
    for (auto it = first; it != last; ++it, ++weight_it) {
        double value = static_cast<double>(*it);
        double weight = static_cast<double>(*weight_it);
 
        if (weight < 0.0) {
            throw std::invalid_argument("statcpp::weighted_mean: negative weight");
        }
 
        sum_weighted += value * weight;
        sum_weights += weight;
    }
 
    if (sum_weights == 0.0) {
        throw std::invalid_argument("statcpp::weighted_mean: sum of weights is zero");
    }
 
    return sum_weighted / sum_weights;
}
 
template <typename Iterator, typename WeightIterator>
[[deprecated("Use weighted_mean(first, last, weight_first, weight_last) overload for range safety")]]
double weighted_mean(Iterator first, Iterator last, WeightIterator weight_first)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_mean: empty range");
    }
 
    double sum_weighted = 0.0;
    double sum_weights = 0.0;
    auto weight_it = weight_first;
 
    for (auto it = first; it != last; ++it, ++weight_it) {
        double value = static_cast<double>(*it);
        double weight = static_cast<double>(*weight_it);
 
        if (weight < 0.0) {
            throw std::invalid_argument("statcpp::weighted_mean: negative weight");
        }
 
        sum_weighted += value * weight;
        sum_weights += weight;
    }
 
    if (sum_weights == 0.0) {
        throw std::invalid_argument("statcpp::weighted_mean: sum of weights is zero");
    }
 
    return sum_weighted / sum_weights;
}
 
template <typename Iterator, typename WeightIterator, typename Projection>
double weighted_mean(Iterator first, Iterator last, WeightIterator weight_first, WeightIterator weight_last, Projection proj)
{
    if (std::distance(first, last) != std::distance(weight_first, weight_last)) {
        throw std::invalid_argument("statcpp::weighted_mean: data and weight ranges must have the same size");
    }
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_mean: empty range");
    }
 
    double sum_weighted = 0.0;
    double sum_weights = 0.0;
    auto weight_it = weight_first;
 
    for (auto it = first; it != last; ++it, ++weight_it) {
        double value = static_cast<double>(std::invoke(proj, *it));
        double weight = static_cast<double>(*weight_it);
 
        if (weight < 0.0) {
            throw std::invalid_argument("statcpp::weighted_mean: negative weight");
        }
 
        sum_weighted += value * weight;
        sum_weights += weight;
    }
 
    if (sum_weights == 0.0) {
        throw std::invalid_argument("statcpp::weighted_mean: sum of weights is zero");
    }
 
    return sum_weighted / sum_weights;
}
 
template <typename Iterator, typename WeightIterator, typename Projection>
[[deprecated("Use weighted_mean(first, last, weight_first, weight_last, proj) overload for range safety")]]
double weighted_mean(Iterator first, Iterator last, WeightIterator weight_first, Projection proj)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_mean: empty range");
    }
 
    double sum_weighted = 0.0;
    double sum_weights = 0.0;
    auto weight_it = weight_first;
 
    for (auto it = first; it != last; ++it, ++weight_it) {
        double value = static_cast<double>(std::invoke(proj, *it));
        double weight = static_cast<double>(*weight_it);
 
        if (weight < 0.0) {
            throw std::invalid_argument("statcpp::weighted_mean: negative weight");
        }
 
        sum_weighted += value * weight;
        sum_weights += weight;
    }
 
    if (sum_weights == 0.0) {
        throw std::invalid_argument("statcpp::weighted_mean: sum of weights is zero");
    }
 
    return sum_weighted / sum_weights;
}
 
// ============================================================================
// Logarithmic Mean
// ============================================================================
 
template <typename T1, typename T2>
double logarithmic_mean(T1 a, T2 b)
{
    double x = static_cast<double>(a);
    double y = static_cast<double>(b);
 
    if (x <= 0.0 || y <= 0.0) {
        throw std::invalid_argument("statcpp::logarithmic_mean: arguments must be positive");
    }
 
    // If values are approximately equal, return x directly.
    // Use relative difference to avoid misclassifying large values as distinct
    // when the relative gap is tiny (e.g., x = 1e15, y = 1e15 + 1e5).
    if (std::abs(x - y) <= 1e-10 * std::max(x, y)) {
        return x;
    }
 
    // Standard logarithmic mean: (y - x) / (ln(y) - ln(x))
    return (y - x) / (std::log(y) - std::log(x));
}
 
// ============================================================================
// Weighted Harmonic Mean
// ============================================================================
 
template <typename Iterator, typename WeightIterator>
double weighted_harmonic_mean(Iterator first, Iterator last, WeightIterator weight_first, WeightIterator weight_last)
{
    if (std::distance(first, last) != std::distance(weight_first, weight_last)) {
        throw std::invalid_argument("statcpp::weighted_harmonic_mean: data and weight ranges must have the same size");
    }
 
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_harmonic_mean: empty range");
    }
 
    double sum_weighted = 0.0;
    double sum_weights = 0.0;
    auto weight_it = weight_first;
 
    for (auto it = first; it != last; ++it, ++weight_it) {
        double value = static_cast<double>(*it);
        double weight = static_cast<double>(*weight_it);
 
        if (weight < 0.0) {
            throw std::invalid_argument("statcpp::weighted_harmonic_mean: negative weight");
        }
 
        if (std::abs(value) < std::numeric_limits<double>::min()) {
            throw std::invalid_argument("statcpp::weighted_harmonic_mean: zero value");
        }
 
        sum_weighted += weight / value;
        sum_weights += weight;
    }
 
    if (sum_weights == 0.0) {
        throw std::invalid_argument("statcpp::weighted_harmonic_mean: sum of weights is zero");
    }
 
    return sum_weights / sum_weighted;
}
 
template <typename Iterator, typename WeightIterator>
[[deprecated("Use weighted_harmonic_mean(first, last, weight_first, weight_last) overload for range safety")]]
double weighted_harmonic_mean(Iterator first, Iterator last, WeightIterator weight_first)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_harmonic_mean: empty range");
    }
 
    double sum_weighted = 0.0;
    double sum_weights = 0.0;
    auto weight_it = weight_first;
 
    for (auto it = first; it != last; ++it, ++weight_it) {
        double value = static_cast<double>(*it);
        double weight = static_cast<double>(*weight_it);
 
        if (weight < 0.0) {
            throw std::invalid_argument("statcpp::weighted_harmonic_mean: negative weight");
        }
 
        if (std::abs(value) < std::numeric_limits<double>::min()) {
            throw std::invalid_argument("statcpp::weighted_harmonic_mean: zero value");
        }
 
        sum_weighted += weight / value;
        sum_weights += weight;
    }
 
    if (sum_weights == 0.0) {
        throw std::invalid_argument("statcpp::weighted_harmonic_mean: sum of weights is zero");
    }
 
    return sum_weights / sum_weighted;
}
 
template <typename Iterator, typename WeightIterator, typename Projection>
double weighted_harmonic_mean(Iterator first, Iterator last, WeightIterator weight_first, WeightIterator weight_last, Projection proj)
{
    if (std::distance(first, last) != std::distance(weight_first, weight_last)) {
        throw std::invalid_argument("statcpp::weighted_harmonic_mean: data and weight ranges must have the same size");
    }
 
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_harmonic_mean: empty range");
    }
 
    double sum_weighted = 0.0;
    double sum_weights = 0.0;
    auto weight_it = weight_first;
 
    for (auto it = first; it != last; ++it, ++weight_it) {
        double value = static_cast<double>(std::invoke(proj, *it));
        double weight = static_cast<double>(*weight_it);
 
        if (weight < 0.0) {
            throw std::invalid_argument("statcpp::weighted_harmonic_mean: negative weight");
        }
 
        if (std::abs(value) < std::numeric_limits<double>::min()) {
            throw std::invalid_argument("statcpp::weighted_harmonic_mean: zero value");
        }
 
        sum_weighted += weight / value;
        sum_weights += weight;
    }
 
    if (sum_weights == 0.0) {
        throw std::invalid_argument("statcpp::weighted_harmonic_mean: sum of weights is zero");
    }
 
    return sum_weights / sum_weighted;
}
 
template <typename Iterator, typename WeightIterator, typename Projection>
[[deprecated("Use weighted_harmonic_mean(first, last, weight_first, weight_last, proj) overload for range safety")]]
double weighted_harmonic_mean(Iterator first, Iterator last, WeightIterator weight_first, Projection proj)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_harmonic_mean: empty range");
    }
 
    double sum_weighted = 0.0;
    double sum_weights = 0.0;
    auto weight_it = weight_first;
 
    for (auto it = first; it != last; ++it, ++weight_it) {
        double value = static_cast<double>(std::invoke(proj, *it));
        double weight = static_cast<double>(*weight_it);
 
        if (weight < 0.0) {
            throw std::invalid_argument("statcpp::weighted_harmonic_mean: negative weight");
        }
 
        if (std::abs(value) < std::numeric_limits<double>::min()) {
            throw std::invalid_argument("statcpp::weighted_harmonic_mean: zero value");
        }
 
        sum_weighted += weight / value;
        sum_weights += weight;
    }
 
    if (sum_weights == 0.0) {
        throw std::invalid_argument("statcpp::weighted_harmonic_mean: sum of weights is zero");
    }
 
    return sum_weights / sum_weighted;
}
 
// ============================================================================
// Argmin / Argmax
// ============================================================================
 
template <typename Iterator>
std::size_t argmin(Iterator first, Iterator last)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::argmin: empty range");
    }
 
    auto min_it = std::min_element(first, last);
    return static_cast<std::size_t>(std::distance(first, min_it));
}
 
template <typename Iterator, typename Projection>
std::size_t argmin(Iterator first, Iterator last, Projection proj)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::argmin: empty range");
    }
 
    auto min_it = std::min_element(first, last,
        [&proj](const auto& a, const auto& b) {
            return std::invoke(proj, a) < std::invoke(proj, b);
        });
 
    return static_cast<std::size_t>(std::distance(first, min_it));
}
 
template <typename Iterator>
std::size_t argmax(Iterator first, Iterator last)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::argmax: empty range");
    }
 
    auto max_it = std::max_element(first, last);
    return static_cast<std::size_t>(std::distance(first, max_it));
}
 
template <typename Iterator, typename Projection>
std::size_t argmax(Iterator first, Iterator last, Projection proj)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::argmax: empty range");
    }
 
    auto max_it = std::max_element(first, last,
        [&proj](const auto& a, const auto& b) {
            return std::invoke(proj, a) < std::invoke(proj, b);
        });
 
    return static_cast<std::size_t>(std::distance(first, max_it));
}
 
} // namespace statcpp