order_statistics.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <functional>
#include <iterator>
#include <type_traits>
#include <limits>
#include <stdexcept>
#include <utility>
#include <vector>
 
namespace statcpp {
 
// ============================================================================
// Result Structures
// ============================================================================
 
struct quartile_result {
    double q1;  
    double q2;  
    double q3;  
};
 
struct five_number_summary_result {
    double min;     
    double q1;      
    double median;  
    double q3;      
    double max;     
};
 
// ============================================================================
// Linear Interpolation Helper (R type=7 / Excel QUARTILE.INC equivalent)
// ============================================================================
 
template <typename Iterator>
double interpolate_at(Iterator first, std::size_t n, double p)
{
    static_assert(
        std::is_base_of_v<
            std::random_access_iterator_tag,
            typename std::iterator_traits<Iterator>::iterator_category>,
        "statcpp::interpolate_at requires random access iterators");
 
    double index = p * static_cast<double>(n - 1);
    auto lo = static_cast<std::size_t>(std::floor(index));
    double frac = index - static_cast<double>(lo);
    if (lo + 1 >= n) {
        return static_cast<double>(*(first + lo));
    }
    return static_cast<double>(*(first + lo)) * (1.0 - frac)
         + static_cast<double>(*(first + lo + 1)) * frac;
}
 
template <typename Iterator, typename Projection>
double interpolate_at(Iterator first, std::size_t n, double p, Projection proj)
{
    static_assert(
        std::is_base_of_v<
            std::random_access_iterator_tag,
            typename std::iterator_traits<Iterator>::iterator_category>,
        "statcpp::interpolate_at requires random access iterators");
 
    double index = p * static_cast<double>(n - 1);
    auto lo = static_cast<std::size_t>(std::floor(index));
    double frac = index - static_cast<double>(lo);
    if (lo + 1 >= n) {
        return static_cast<double>(std::invoke(proj, *(first + lo)));
    }
    return static_cast<double>(std::invoke(proj, *(first + lo))) * (1.0 - frac)
         + static_cast<double>(std::invoke(proj, *(first + lo + 1))) * frac;
}
 
// ============================================================================
// Minimum
// ============================================================================
 
template <typename Iterator>
auto minimum(Iterator first, Iterator last)
{
    if (first == last) {
        throw std::invalid_argument("statcpp::minimum: empty range");
    }
    return *std::min_element(first, last);
}
 
template <typename Iterator, typename Projection>
auto minimum(Iterator first, Iterator last, Projection proj)
{
    if (first == last) {
        throw std::invalid_argument("statcpp::minimum: empty range");
    }
    auto min_val = std::invoke(proj, *first);
    for (auto it = std::next(first); it != last; ++it) {
        auto val = std::invoke(proj, *it);
        if (val < min_val) {
            min_val = val;
        }
    }
    return min_val;
}
 
// ============================================================================
// Maximum
// ============================================================================
 
template <typename Iterator>
auto maximum(Iterator first, Iterator last)
{
    if (first == last) {
        throw std::invalid_argument("statcpp::maximum: empty range");
    }
    return *std::max_element(first, last);
}
 
template <typename Iterator, typename Projection>
auto maximum(Iterator first, Iterator last, Projection proj)
{
    if (first == last) {
        throw std::invalid_argument("statcpp::maximum: empty range");
    }
    auto max_val = std::invoke(proj, *first);
    for (auto it = std::next(first); it != last; ++it) {
        auto val = std::invoke(proj, *it);
        if (val > max_val) {
            max_val = val;
        }
    }
    return max_val;
}
 
// ============================================================================
// Quartiles (Q1, Q2, Q3)
// ============================================================================
 
template <typename Iterator>
quartile_result quartiles(Iterator first, Iterator last)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::quartiles: empty range");
    }
    return {
        interpolate_at(first, n, 0.25),
        interpolate_at(first, n, 0.50),
        interpolate_at(first, n, 0.75)
    };
}
 
template <typename Iterator, typename Projection>
quartile_result quartiles(Iterator first, Iterator last, Projection proj)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::quartiles: empty range");
    }
    return {
        interpolate_at(first, n, 0.25, proj),
        interpolate_at(first, n, 0.50, proj),
        interpolate_at(first, n, 0.75, proj)
    };
}
 
// ============================================================================
// Percentile
// ============================================================================
 
template <typename Iterator>
double percentile(Iterator first, Iterator last, double p)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::percentile: empty range");
    }
    if (p < 0.0 || p > 1.0) {
        throw std::invalid_argument("statcpp::percentile: p must be in [0, 1]");
    }
    return interpolate_at(first, n, p);
}
 
template <typename Iterator, typename Projection>
double percentile(Iterator first, Iterator last, double p, Projection proj)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::percentile: empty range");
    }
    if (p < 0.0 || p > 1.0) {
        throw std::invalid_argument("statcpp::percentile: p must be in [0, 1]");
    }
    return interpolate_at(first, n, p, proj);
}
 
// ============================================================================
// Five-Number Summary
// ============================================================================
 
template <typename Iterator>
five_number_summary_result five_number_summary(Iterator first, Iterator last)
{
    static_assert(
        std::is_base_of_v<
            std::random_access_iterator_tag,
            typename std::iterator_traits<Iterator>::iterator_category>,
        "statcpp::five_number_summary requires random access iterators");
 
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::five_number_summary: empty range");
    }
    return {
        static_cast<double>(*first),
        interpolate_at(first, n, 0.25),
        interpolate_at(first, n, 0.50),
        interpolate_at(first, n, 0.75),
        static_cast<double>(*(first + (n - 1)))
    };
}
 
template <typename Iterator, typename Projection>
five_number_summary_result five_number_summary(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::five_number_summary requires random access iterators");
 
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::five_number_summary: empty range");
    }
    return {
        static_cast<double>(std::invoke(proj, *first)),
        interpolate_at(first, n, 0.25, proj),
        interpolate_at(first, n, 0.50, proj),
        interpolate_at(first, n, 0.75, proj),
        static_cast<double>(std::invoke(proj, *(first + (n - 1))))
    };
}
 
// ============================================================================
// Weighted Median
// ============================================================================
 
template <typename Iterator, typename WeightIterator>
double weighted_median(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_median: data and weight ranges must have the same size");
    }
 
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_median: empty range");
    }
 
    // Create value-weight pairs
    std::vector<std::pair<double, double>> pairs;
    pairs.reserve(n);
    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_median: negative weight");
        }
        pairs.emplace_back(value, weight);
    }
 
    // Sort by value
    std::sort(pairs.begin(), pairs.end(),
              [](const auto& a, const auto& b) { return a.first < b.first; });
 
    // Sum of weights
    double total_weight = 0.0;
    for (const auto& p : pairs) {
        total_weight += p.second;
    }
 
    if (total_weight == 0.0) {
        throw std::invalid_argument("statcpp::weighted_median: sum of weights is zero");
    }
 
    // Calculate cumulative weight to find median.
    // Use a relative tolerance when testing whether cumulative == half_weight.
    // Plain == on a floating-point accumulator is unreliable because sequential
    // additions may produce a value within rounding error of half_weight without
    // being bit-identical.
    double cumulative = 0.0;
    double half_weight = total_weight / 2.0;
    const double tol = std::numeric_limits<double>::epsilon() * half_weight;
 
    for (std::size_t i = 0; i < pairs.size(); ++i) {
        cumulative += pairs[i].second;
        if (cumulative >= half_weight) {
            // If cumulative weight is (approximately) exactly half, take average with next positive-weight value
            if (std::abs(cumulative - half_weight) <= tol) {
                for (std::size_t j = i + 1; j < pairs.size(); ++j) {
                    if (pairs[j].second > 0.0) {
                        return (pairs[i].first + pairs[j].first) / 2.0;
                    }
                }
            }
            return pairs[i].first;
        }
    }
 
    return pairs.back().first;
}
 
template <typename Iterator, typename WeightIterator>
[[deprecated("Use weighted_median(first, last, weight_first, weight_last) overload for range safety")]]
double weighted_median(Iterator first, Iterator last, WeightIterator weight_first)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_median: empty range");
    }
 
    // Create value-weight pairs
    std::vector<std::pair<double, double>> pairs;
    pairs.reserve(n);
    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_median: negative weight");
        }
        pairs.emplace_back(value, weight);
    }
 
    // Sort by value
    std::sort(pairs.begin(), pairs.end(),
              [](const auto& a, const auto& b) { return a.first < b.first; });
 
    // Sum of weights
    double total_weight = 0.0;
    for (const auto& p : pairs) {
        total_weight += p.second;
    }
 
    if (total_weight == 0.0) {
        throw std::invalid_argument("statcpp::weighted_median: sum of weights is zero");
    }
 
    // Calculate cumulative weight to find median.
    // Use a relative tolerance when testing whether cumulative == half_weight.
    // Plain == on a floating-point accumulator is unreliable because sequential
    // additions may produce a value within rounding error of half_weight without
    // being bit-identical.
    double cumulative = 0.0;
    double half_weight = total_weight / 2.0;
    const double tol = std::numeric_limits<double>::epsilon() * half_weight;
 
    for (std::size_t i = 0; i < pairs.size(); ++i) {
        cumulative += pairs[i].second;
        if (cumulative >= half_weight) {
            // If cumulative weight is (approximately) exactly half, take average with next positive-weight value
            if (std::abs(cumulative - half_weight) <= tol) {
                for (std::size_t j = i + 1; j < pairs.size(); ++j) {
                    if (pairs[j].second > 0.0) {
                        return (pairs[i].first + pairs[j].first) / 2.0;
                    }
                }
            }
            return pairs[i].first;
        }
    }
 
    return pairs.back().first;
}
 
template <typename Iterator, typename WeightIterator, typename Projection>
double weighted_median(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_median: data and weight ranges must have the same size");
    }
 
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_median: empty range");
    }
 
    // Create value-weight pairs
    std::vector<std::pair<double, double>> pairs;
    pairs.reserve(n);
    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_median: negative weight");
        }
        pairs.emplace_back(value, weight);
    }
 
    // Sort by value
    std::sort(pairs.begin(), pairs.end(),
              [](const auto& a, const auto& b) { return a.first < b.first; });
 
    // Sum of weights
    double total_weight = 0.0;
    for (const auto& p : pairs) {
        total_weight += p.second;
    }
 
    if (total_weight == 0.0) {
        throw std::invalid_argument("statcpp::weighted_median: sum of weights is zero");
    }
 
    // Calculate cumulative weight to find median.
    // Use a relative tolerance when testing whether cumulative == half_weight.
    // Plain == on a floating-point accumulator is unreliable because sequential
    // additions may produce a value within rounding error of half_weight without
    // being bit-identical.
    double cumulative = 0.0;
    double half_weight = total_weight / 2.0;
    const double tol = std::numeric_limits<double>::epsilon() * half_weight;
 
    for (std::size_t i = 0; i < pairs.size(); ++i) {
        cumulative += pairs[i].second;
        if (cumulative >= half_weight) {
            // If cumulative weight is (approximately) exactly half, take average with next positive-weight value
            if (std::abs(cumulative - half_weight) <= tol) {
                for (std::size_t j = i + 1; j < pairs.size(); ++j) {
                    if (pairs[j].second > 0.0) {
                        return (pairs[i].first + pairs[j].first) / 2.0;
                    }
                }
            }
            return pairs[i].first;
        }
    }
 
    return pairs.back().first;
}
 
template <typename Iterator, typename WeightIterator, typename Projection>
[[deprecated("Use weighted_median(first, last, weight_first, weight_last, proj) overload for range safety")]]
double weighted_median(Iterator first, Iterator last, WeightIterator weight_first, Projection proj)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_median: empty range");
    }
 
    // Create value-weight pairs
    std::vector<std::pair<double, double>> pairs;
    pairs.reserve(n);
    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_median: negative weight");
        }
        pairs.emplace_back(value, weight);
    }
 
    // Sort by value
    std::sort(pairs.begin(), pairs.end(),
              [](const auto& a, const auto& b) { return a.first < b.first; });
 
    // Sum of weights
    double total_weight = 0.0;
    for (const auto& p : pairs) {
        total_weight += p.second;
    }
 
    if (total_weight == 0.0) {
        throw std::invalid_argument("statcpp::weighted_median: sum of weights is zero");
    }
 
    // Calculate cumulative weight to find median.
    // Use a relative tolerance when testing whether cumulative == half_weight.
    // Plain == on a floating-point accumulator is unreliable because sequential
    // additions may produce a value within rounding error of half_weight without
    // being bit-identical.
    double cumulative = 0.0;
    double half_weight = total_weight / 2.0;
    const double tol = std::numeric_limits<double>::epsilon() * half_weight;
 
    for (std::size_t i = 0; i < pairs.size(); ++i) {
        cumulative += pairs[i].second;
        if (cumulative >= half_weight) {
            // If cumulative weight is (approximately) exactly half, take average with next positive-weight value
            if (std::abs(cumulative - half_weight) <= tol) {
                for (std::size_t j = i + 1; j < pairs.size(); ++j) {
                    if (pairs[j].second > 0.0) {
                        return (pairs[i].first + pairs[j].first) / 2.0;
                    }
                }
            }
            return pairs[i].first;
        }
    }
 
    return pairs.back().first;
}
 
// ============================================================================
// Weighted Percentile
// ============================================================================
 
template <typename Iterator, typename WeightIterator>
double weighted_percentile(Iterator first, Iterator last, WeightIterator weight_first, WeightIterator weight_last, double p)
{
    if (std::distance(first, last) != std::distance(weight_first, weight_last)) {
        throw std::invalid_argument("statcpp::weighted_percentile: data and weight ranges must have the same size");
    }
 
    if (!(0.0 <= p && p <= 1.0)) {
        throw std::invalid_argument("statcpp::weighted_percentile: p must be in [0, 1]");
    }
 
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_percentile: empty range");
    }
 
    // Create value-weight pairs
    std::vector<std::pair<double, double>> pairs;
    pairs.reserve(n);
    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_percentile: negative weight");
        }
        pairs.emplace_back(value, weight);
    }
 
    // Sort by value
    std::sort(pairs.begin(), pairs.end(),
              [](const auto& a, const auto& b) { return a.first < b.first; });
 
    // Sum of weights
    double total_weight = 0.0;
    for (const auto& pair : pairs) {
        total_weight += pair.second;
    }
 
    if (total_weight == 0.0) {
        throw std::invalid_argument("statcpp::weighted_percentile: sum of weights is zero");
    }
 
    // Explicit endpoint handling to avoid relying on floating-point loop
    if (p <= 0.0) return pairs.front().first;
    if (p >= 1.0) return pairs.back().first;
 
    // Calculate cumulative weight to find target percentile
    // For p=0.5, find position at 50% of total weight
    double target = p * total_weight;
    const double tol = std::numeric_limits<double>::epsilon() * total_weight;
    double cumulative = 0.0;
 
    for (std::size_t i = 0; i < pairs.size(); ++i) {
        cumulative += pairs[i].second;
        if (cumulative >= target) {
            // If cumulative weight is approximately at target, take average with next value
            if (std::abs(cumulative - target) <= tol && i + 1 < pairs.size()) {
                return (pairs[i].first + pairs[i + 1].first) / 2.0;
            }
            return pairs[i].first;
        }
    }
 
    return pairs.back().first;
}
 
template <typename Iterator, typename WeightIterator>
[[deprecated("Use weighted_percentile(first, last, weight_first, weight_last, p) overload for range safety")]]
double weighted_percentile(Iterator first, Iterator last, WeightIterator weight_first, double p)
{
    if (!(0.0 <= p && p <= 1.0)) {
        throw std::invalid_argument("statcpp::weighted_percentile: p must be in [0, 1]");
    }
 
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_percentile: empty range");
    }
 
    // Create value-weight pairs
    std::vector<std::pair<double, double>> pairs;
    pairs.reserve(n);
    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_percentile: negative weight");
        }
        pairs.emplace_back(value, weight);
    }
 
    // Sort by value
    std::sort(pairs.begin(), pairs.end(),
              [](const auto& a, const auto& b) { return a.first < b.first; });
 
    // Sum of weights
    double total_weight = 0.0;
    for (const auto& pair : pairs) {
        total_weight += pair.second;
    }
 
    if (total_weight == 0.0) {
        throw std::invalid_argument("statcpp::weighted_percentile: sum of weights is zero");
    }
 
    // Explicit endpoint handling to avoid relying on floating-point loop
    if (p <= 0.0) return pairs.front().first;
    if (p >= 1.0) return pairs.back().first;
 
    // Calculate cumulative weight to find target percentile
    // For p=0.5, find position at 50% of total weight
    double target = p * total_weight;
    const double tol = std::numeric_limits<double>::epsilon() * total_weight;
    double cumulative = 0.0;
 
    for (std::size_t i = 0; i < pairs.size(); ++i) {
        cumulative += pairs[i].second;
        if (cumulative >= target) {
            // If cumulative weight is approximately at target, take average with next value
            if (std::abs(cumulative - target) <= tol && i + 1 < pairs.size()) {
                return (pairs[i].first + pairs[i + 1].first) / 2.0;
            }
            return pairs[i].first;
        }
    }
 
    return pairs.back().first;
}
 
template <typename Iterator, typename WeightIterator, typename Projection>
double weighted_percentile(Iterator first, Iterator last, WeightIterator weight_first, WeightIterator weight_last, double p, Projection proj)
{
    if (std::distance(first, last) != std::distance(weight_first, weight_last)) {
        throw std::invalid_argument("statcpp::weighted_percentile: data and weight ranges must have the same size");
    }
 
    if (!(0.0 <= p && p <= 1.0)) {
        throw std::invalid_argument("statcpp::weighted_percentile: p must be in [0, 1]");
    }
 
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_percentile: empty range");
    }
 
    // Create value-weight pairs
    std::vector<std::pair<double, double>> pairs;
    pairs.reserve(n);
    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_percentile: negative weight");
        }
        pairs.emplace_back(value, weight);
    }
 
    // Sort by value
    std::sort(pairs.begin(), pairs.end(),
              [](const auto& a, const auto& b) { return a.first < b.first; });
 
    // Sum of weights
    double total_weight = 0.0;
    for (const auto& pair : pairs) {
        total_weight += pair.second;
    }
 
    if (total_weight == 0.0) {
        throw std::invalid_argument("statcpp::weighted_percentile: sum of weights is zero");
    }
 
    // Explicit endpoint handling to avoid relying on floating-point loop
    if (p <= 0.0) return pairs.front().first;
    if (p >= 1.0) return pairs.back().first;
 
    // Calculate cumulative weight to find target percentile
    // For p=0.5, find position at 50% of total weight
    double target = p * total_weight;
    const double tol = std::numeric_limits<double>::epsilon() * total_weight;
    double cumulative = 0.0;
 
    for (std::size_t i = 0; i < pairs.size(); ++i) {
        cumulative += pairs[i].second;
        if (cumulative >= target) {
            // If cumulative weight is approximately at target, take average with next value
            if (std::abs(cumulative - target) <= tol && i + 1 < pairs.size()) {
                return (pairs[i].first + pairs[i + 1].first) / 2.0;
            }
            return pairs[i].first;
        }
    }
 
    return pairs.back().first;
}
 
template <typename Iterator, typename WeightIterator, typename Projection>
[[deprecated("Use weighted_percentile(first, last, weight_first, weight_last, p, proj) overload for range safety")]]
double weighted_percentile(Iterator first, Iterator last, WeightIterator weight_first, double p, Projection proj)
{
    if (!(0.0 <= p && p <= 1.0)) {
        throw std::invalid_argument("statcpp::weighted_percentile: p must be in [0, 1]");
    }
 
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_percentile: empty range");
    }
 
    // Create value-weight pairs
    std::vector<std::pair<double, double>> pairs;
    pairs.reserve(n);
    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_percentile: negative weight");
        }
        pairs.emplace_back(value, weight);
    }
 
    // Sort by value
    std::sort(pairs.begin(), pairs.end(),
              [](const auto& a, const auto& b) { return a.first < b.first; });
 
    // Sum of weights
    double total_weight = 0.0;
    for (const auto& pair : pairs) {
        total_weight += pair.second;
    }
 
    if (total_weight == 0.0) {
        throw std::invalid_argument("statcpp::weighted_percentile: sum of weights is zero");
    }
 
    // Explicit endpoint handling to avoid relying on floating-point loop
    if (p <= 0.0) return pairs.front().first;
    if (p >= 1.0) return pairs.back().first;
 
    // Calculate cumulative weight to find target percentile
    // For p=0.5, find position at 50% of total weight
    double target = p * total_weight;
    const double tol = std::numeric_limits<double>::epsilon() * total_weight;
    double cumulative = 0.0;
 
    for (std::size_t i = 0; i < pairs.size(); ++i) {
        cumulative += pairs[i].second;
        if (cumulative >= target) {
            // If cumulative weight is approximately at target, take average with next value
            if (std::abs(cumulative - target) <= tol && i + 1 < pairs.size()) {
                return (pairs[i].first + pairs[i + 1].first) / 2.0;
            }
            return pairs[i].first;
        }
    }
 
    return pairs.back().first;
}
 
} // namespace statcpp