dispersion_spread.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <functional>
#include <iterator>
#include <stdexcept>
#include <type_traits>
#include <vector>
 
#include "statcpp/basic_statistics.hpp"
#include "statcpp/order_statistics.hpp"
 
namespace statcpp {
 
// ============================================================================
// Range
// ============================================================================
 
template <typename Iterator>
double range(Iterator first, Iterator last)
{
    if (first == last) {
        throw std::invalid_argument("statcpp::range: empty range");
    }
    auto [min_it, max_it] = std::minmax_element(first, last);
    return static_cast<double>(*max_it) - static_cast<double>(*min_it);
}
 
template <typename Iterator, typename Projection>
double range(Iterator first, Iterator last, Projection proj)
{
    if (first == last) {
        throw std::invalid_argument("statcpp::range: empty range");
    }
    auto it = first;
    double min_val = static_cast<double>(std::invoke(proj, *it));
    double max_val = min_val;
    ++it;
    for (; it != last; ++it) {
        double val = static_cast<double>(std::invoke(proj, *it));
        if (val < min_val) {
            min_val = val;
        }
        if (val > max_val) {
            max_val = val;
        }
    }
    return max_val - min_val;
}
 
// ============================================================================
// Variance with ddof
// ============================================================================
 
template <typename Iterator>
double var(Iterator first, Iterator last, std::size_t ddof = 0)
{
    if (ddof > 1) {
        throw std::invalid_argument("statcpp::var: ddof must be 0 or 1");
    }
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::var: empty range");
    }
    if (ddof == 1 && n < 2) {
        throw std::invalid_argument("statcpp::var: need at least 2 elements for ddof=1");
    }
    double m = statcpp::mean(first, last);
    double sum_sq = 0.0;
    for (auto it = first; it != last; ++it) {
        double diff = static_cast<double>(*it) - m;
        sum_sq += diff * diff;
    }
    return sum_sq / static_cast<double>(n - ddof);
}
 
template <typename Iterator>
double var(Iterator first, Iterator last, double precomputed_mean, std::size_t ddof)
{
    if (ddof > 1) {
        throw std::invalid_argument("statcpp::var: ddof must be 0 or 1");
    }
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::var: empty range");
    }
    if (ddof == 1 && n < 2) {
        throw std::invalid_argument("statcpp::var: need at least 2 elements for ddof=1");
    }
    double sum_sq = 0.0;
    for (auto it = first; it != last; ++it) {
        double diff = static_cast<double>(*it) - precomputed_mean;
        sum_sq += diff * diff;
    }
    return sum_sq / static_cast<double>(n - ddof);
}
 
template <typename Iterator, typename Projection,
          typename = std::enable_if_t<
              std::is_invocable_v<Projection,
                  typename std::iterator_traits<Iterator>::value_type>>>
double var(Iterator first, Iterator last, Projection proj, std::size_t ddof = 0)
{
    if (ddof > 1) {
        throw std::invalid_argument("statcpp::var: ddof must be 0 or 1");
    }
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::var: empty range");
    }
    if (ddof == 1 && n < 2) {
        throw std::invalid_argument("statcpp::var: need at least 2 elements for ddof=1");
    }
    double m = statcpp::mean(first, last, proj);
    double sum_sq = 0.0;
    for (auto it = first; it != last; ++it) {
        double diff = static_cast<double>(std::invoke(proj, *it)) - m;
        sum_sq += diff * diff;
    }
    return sum_sq / static_cast<double>(n - ddof);
}
 
template <typename Iterator, typename Projection>
double var(Iterator first, Iterator last, Projection proj, double precomputed_mean, std::size_t ddof)
{
    if (ddof > 1) {
        throw std::invalid_argument("statcpp::var: ddof must be 0 or 1");
    }
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::var: empty range");
    }
    if (ddof == 1 && n < 2) {
        throw std::invalid_argument("statcpp::var: need at least 2 elements for ddof=1");
    }
    double sum_sq = 0.0;
    for (auto it = first; it != last; ++it) {
        double diff = static_cast<double>(std::invoke(proj, *it)) - precomputed_mean;
        sum_sq += diff * diff;
    }
    return sum_sq / static_cast<double>(n - ddof);
}
 
// ============================================================================
// Population Variance (alias for ddof=0)
// ============================================================================
 
template <typename Iterator>
double population_variance(Iterator first, Iterator last)
{
    return var(first, last, 0);
}
 
template <typename Iterator>
double population_variance(Iterator first, Iterator last, double precomputed_mean)
{
    return var(first, last, precomputed_mean, 0);
}
 
template <typename Iterator, typename Projection,
          typename = std::enable_if_t<
              std::is_invocable_v<Projection,
                  typename std::iterator_traits<Iterator>::value_type>>>
double population_variance(Iterator first, Iterator last, Projection proj)
{
    return var(first, last, proj, 0);
}
 
template <typename Iterator, typename Projection>
double population_variance(Iterator first, Iterator last, Projection proj, double precomputed_mean)
{
    return var(first, last, proj, precomputed_mean, 0);
}
 
// ============================================================================
// Sample Variance (alias for ddof=1)
// ============================================================================
 
template <typename Iterator>
double sample_variance(Iterator first, Iterator last)
{
    return var(first, last, 1);
}
 
template <typename Iterator>
double sample_variance(Iterator first, Iterator last, double precomputed_mean)
{
    return var(first, last, precomputed_mean, 1);
}
 
template <typename Iterator, typename Projection,
          typename = std::enable_if_t<
              std::is_invocable_v<Projection,
                  typename std::iterator_traits<Iterator>::value_type>>>
double sample_variance(Iterator first, Iterator last, Projection proj)
{
    return var(first, last, proj, 1);
}
 
template <typename Iterator, typename Projection>
double sample_variance(Iterator first, Iterator last, Projection proj, double precomputed_mean)
{
    return var(first, last, proj, precomputed_mean, 1);
}
 
// ============================================================================
// Variance (= Sample Variance)
// ============================================================================
 
template <typename Iterator>
double variance(Iterator first, Iterator last)
{
    return sample_variance(first, last);
}
 
template <typename Iterator>
double variance(Iterator first, Iterator last, double precomputed_mean)
{
    return sample_variance(first, last, precomputed_mean);
}
 
template <typename Iterator, typename Projection,
          typename = std::enable_if_t<
              std::is_invocable_v<Projection,
                  typename std::iterator_traits<Iterator>::value_type>>>
double variance(Iterator first, Iterator last, Projection proj)
{
    return sample_variance(first, last, proj);
}
 
template <typename Iterator, typename Projection>
double variance(Iterator first, Iterator last, Projection proj, double precomputed_mean)
{
    return sample_variance(first, last, proj, precomputed_mean);
}
 
// ============================================================================
// Standard Deviation with ddof (NumPy-style)
// ============================================================================
 
template <typename Iterator>
double stdev(Iterator first, Iterator last, std::size_t ddof = 0)
{
    return std::sqrt(var(first, last, ddof));
}
 
template <typename Iterator>
double stdev(Iterator first, Iterator last, double precomputed_mean, std::size_t ddof)
{
    return std::sqrt(var(first, last, precomputed_mean, ddof));
}
 
template <typename Iterator, typename Projection,
          typename = std::enable_if_t<
              std::is_invocable_v<Projection,
                  typename std::iterator_traits<Iterator>::value_type>>>
double stdev(Iterator first, Iterator last, Projection proj, std::size_t ddof = 0)
{
    return std::sqrt(var(first, last, proj, ddof));
}
 
template <typename Iterator, typename Projection>
double stdev(Iterator first, Iterator last, Projection proj, double precomputed_mean, std::size_t ddof)
{
    return std::sqrt(var(first, last, proj, precomputed_mean, ddof));
}
 
// ============================================================================
// Population Standard Deviation (alias for ddof=0)
// ============================================================================
 
template <typename Iterator>
double population_stddev(Iterator first, Iterator last)
{
    return stdev(first, last, 0);
}
 
template <typename Iterator>
double population_stddev(Iterator first, Iterator last, double precomputed_mean)
{
    return stdev(first, last, precomputed_mean, 0);
}
 
template <typename Iterator, typename Projection,
          typename = std::enable_if_t<
              std::is_invocable_v<Projection,
                  typename std::iterator_traits<Iterator>::value_type>>>
double population_stddev(Iterator first, Iterator last, Projection proj)
{
    return stdev(first, last, proj, 0);
}
 
template <typename Iterator, typename Projection>
double population_stddev(Iterator first, Iterator last, Projection proj, double precomputed_mean)
{
    return stdev(first, last, proj, precomputed_mean, 0);
}
 
// ============================================================================
// Sample Standard Deviation (alias for ddof=1)
// ============================================================================
 
template <typename Iterator>
double sample_stddev(Iterator first, Iterator last)
{
    return stdev(first, last, 1);
}
 
template <typename Iterator>
double sample_stddev(Iterator first, Iterator last, double precomputed_mean)
{
    return stdev(first, last, precomputed_mean, 1);
}
 
template <typename Iterator, typename Projection,
          typename = std::enable_if_t<
              std::is_invocable_v<Projection,
                  typename std::iterator_traits<Iterator>::value_type>>>
double sample_stddev(Iterator first, Iterator last, Projection proj)
{
    return stdev(first, last, proj, 1);
}
 
template <typename Iterator, typename Projection>
double sample_stddev(Iterator first, Iterator last, Projection proj, double precomputed_mean)
{
    return stdev(first, last, proj, precomputed_mean, 1);
}
 
// ============================================================================
// Standard Deviation (= Sample Standard Deviation)
// ============================================================================
 
template <typename Iterator>
double stddev(Iterator first, Iterator last)
{
    return sample_stddev(first, last);
}
 
template <typename Iterator>
double stddev(Iterator first, Iterator last, double precomputed_mean)
{
    return sample_stddev(first, last, precomputed_mean);
}
 
template <typename Iterator, typename Projection,
          typename = std::enable_if_t<
              std::is_invocable_v<Projection,
                  typename std::iterator_traits<Iterator>::value_type>>>
double stddev(Iterator first, Iterator last, Projection proj)
{
    return sample_stddev(first, last, proj);
}
 
template <typename Iterator, typename Projection>
double stddev(Iterator first, Iterator last, Projection proj, double precomputed_mean)
{
    return sample_stddev(first, last, proj, precomputed_mean);
}
 
// ============================================================================
// Coefficient of Variation
// ============================================================================
 
template <typename Iterator>
double coefficient_of_variation(Iterator first, Iterator last)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::coefficient_of_variation: empty range");
    }
    if (n < 2) {
        throw std::invalid_argument("statcpp::coefficient_of_variation: need at least 2 elements");
    }
    double m = statcpp::mean(first, last);
    if (m == 0.0) {
        throw std::invalid_argument("statcpp::coefficient_of_variation: mean is zero");
    }
    return sample_stddev(first, last, m) / std::abs(m);
}
 
template <typename Iterator>
double coefficient_of_variation(Iterator first, Iterator last, double precomputed_mean)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::coefficient_of_variation: empty range");
    }
    if (n < 2) {
        throw std::invalid_argument("statcpp::coefficient_of_variation: need at least 2 elements");
    }
    if (precomputed_mean == 0.0) {
        throw std::invalid_argument("statcpp::coefficient_of_variation: mean is zero");
    }
    return sample_stddev(first, last, precomputed_mean) / std::abs(precomputed_mean);
}
 
template <typename Iterator, typename Projection,
          typename = std::enable_if_t<
              std::is_invocable_v<Projection,
                  typename std::iterator_traits<Iterator>::value_type>>>
double coefficient_of_variation(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::coefficient_of_variation: empty range");
    }
    if (n < 2) {
        throw std::invalid_argument("statcpp::coefficient_of_variation: need at least 2 elements");
    }
    double m = statcpp::mean(first, last, proj);
    if (m == 0.0) {
        throw std::invalid_argument("statcpp::coefficient_of_variation: mean is zero");
    }
    return sample_stddev(first, last, proj, m) / std::abs(m);
}
 
template <typename Iterator, typename Projection>
double coefficient_of_variation(Iterator first, Iterator last, Projection proj, double precomputed_mean)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::coefficient_of_variation: empty range");
    }
    if (n < 2) {
        throw std::invalid_argument("statcpp::coefficient_of_variation: need at least 2 elements");
    }
    if (precomputed_mean == 0.0) {
        throw std::invalid_argument("statcpp::coefficient_of_variation: mean is zero");
    }
    return sample_stddev(first, last, proj, precomputed_mean) / std::abs(precomputed_mean);
}
 
// ============================================================================
// Interquartile Range (IQR)
// ============================================================================
 
template <typename Iterator>
double iqr(Iterator first, Iterator last)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::iqr: empty range");
    }
    double q1 = statcpp::interpolate_at(first, n, 0.25);
    double q3 = statcpp::interpolate_at(first, n, 0.75);
    return q3 - q1;
}
 
template <typename Iterator, typename Projection>
double iqr(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::iqr: empty range");
    }
    double q1 = statcpp::interpolate_at(first, n, 0.25, proj);
    double q3 = statcpp::interpolate_at(first, n, 0.75, proj);
    return q3 - q1;
}
 
// ============================================================================
// Mean Absolute Deviation
// ============================================================================
 
template <typename Iterator>
double mean_absolute_deviation(Iterator first, Iterator last)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::mean_absolute_deviation: empty range");
    }
    double m = statcpp::mean(first, last);
    double sum_abs = 0.0;
    for (auto it = first; it != last; ++it) {
        sum_abs += std::abs(static_cast<double>(*it) - m);
    }
    return sum_abs / static_cast<double>(n);
}
 
template <typename Iterator>
double mean_absolute_deviation(Iterator first, Iterator last, double precomputed_mean)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::mean_absolute_deviation: empty range");
    }
    double sum_abs = 0.0;
    for (auto it = first; it != last; ++it) {
        sum_abs += std::abs(static_cast<double>(*it) - precomputed_mean);
    }
    return sum_abs / static_cast<double>(n);
}
 
template <typename Iterator, typename Projection,
          typename = std::enable_if_t<
              std::is_invocable_v<Projection,
                  typename std::iterator_traits<Iterator>::value_type>>>
double mean_absolute_deviation(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::mean_absolute_deviation: empty range");
    }
    double m = statcpp::mean(first, last, proj);
    double sum_abs = 0.0;
    for (auto it = first; it != last; ++it) {
        sum_abs += std::abs(static_cast<double>(std::invoke(proj, *it)) - m);
    }
    return sum_abs / static_cast<double>(n);
}
 
template <typename Iterator, typename Projection>
double mean_absolute_deviation(Iterator first, Iterator last, Projection proj, double precomputed_mean)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::mean_absolute_deviation: empty range");
    }
    double sum_abs = 0.0;
    for (auto it = first; it != last; ++it) {
        sum_abs += std::abs(static_cast<double>(std::invoke(proj, *it)) - precomputed_mean);
    }
    return sum_abs / static_cast<double>(n);
}
 
// ============================================================================
// Weighted Variance
// ============================================================================
 
template <typename Iterator, typename WeightIterator>
double weighted_variance(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_variance: data and weight ranges differ in length");
    }
 
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_variance: empty range");
    }
 
    // Compute weighted mean
    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_variance: negative weight");
        }
 
        sum_weighted += value * weight;
        sum_weights += weight;
    }
 
    if (sum_weights == 0.0) {
        throw std::invalid_argument("statcpp::weighted_variance: sum of weights is zero");
    }
 
    double mean = sum_weighted / sum_weights;
 
    // Compute weighted variance (with Bessel's correction)
    double sum_squared_dev = 0.0;
    double sum_weights_squared = 0.0;
    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);
        double dev = value - mean;
 
        sum_squared_dev += weight * dev * dev;
        sum_weights_squared += weight * weight;
    }
 
    // Bessel's correction: V = Σw(x-μ)² / ((Σw)² - Σ(w²)) * Σw
    // Simplified: V = Σw(x-μ)² / (Σw - Σw²/Σw)
    double correction = sum_weights - (sum_weights_squared / sum_weights);
 
    if (correction <= 0.0) {
        throw std::invalid_argument("statcpp::weighted_variance: insufficient effective sample size");
    }
 
    return sum_squared_dev / correction;
}
 
template <typename Iterator, typename WeightIterator>
[[deprecated("Use weighted_variance(first, last, weight_first, weight_last) overload for range safety")]]
double weighted_variance(Iterator first, Iterator last, WeightIterator weight_first)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_variance: empty range");
    }
 
    // Compute weighted mean
    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_variance: negative weight");
        }
 
        sum_weighted += value * weight;
        sum_weights += weight;
    }
 
    if (sum_weights == 0.0) {
        throw std::invalid_argument("statcpp::weighted_variance: sum of weights is zero");
    }
 
    double mean = sum_weighted / sum_weights;
 
    // Compute weighted variance (with Bessel's correction)
    double sum_squared_dev = 0.0;
    double sum_weights_squared = 0.0;
    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);
        double dev = value - mean;
 
        sum_squared_dev += weight * dev * dev;
        sum_weights_squared += weight * weight;
    }
 
    // Bessel's correction: V = Σw(x-μ)² / ((Σw)² - Σ(w²)) * Σw
    // Simplified: V = Σw(x-μ)² / (Σw - Σw²/Σw)
    double correction = sum_weights - (sum_weights_squared / sum_weights);
 
    if (correction <= 0.0) {
        throw std::invalid_argument("statcpp::weighted_variance: insufficient effective sample size");
    }
 
    return sum_squared_dev / correction;
}
 
template <typename Iterator, typename WeightIterator, typename Projection>
double weighted_variance(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_variance: data and weight ranges differ in length");
    }
 
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_variance: empty range");
    }
 
    // Compute weighted mean
    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_variance: negative weight");
        }
 
        sum_weighted += value * weight;
        sum_weights += weight;
    }
 
    if (sum_weights == 0.0) {
        throw std::invalid_argument("statcpp::weighted_variance: sum of weights is zero");
    }
 
    double mean = sum_weighted / sum_weights;
 
    // Compute weighted variance
    double sum_squared_dev = 0.0;
    double sum_weights_squared = 0.0;
    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);
        double dev = value - mean;
 
        sum_squared_dev += weight * dev * dev;
        sum_weights_squared += weight * weight;
    }
 
    double correction = sum_weights - (sum_weights_squared / sum_weights);
 
    if (correction <= 0.0) {
        throw std::invalid_argument("statcpp::weighted_variance: insufficient effective sample size");
    }
 
    return sum_squared_dev / correction;
}
 
template <typename Iterator, typename WeightIterator, typename Projection>
[[deprecated("Use weighted_variance(first, last, weight_first, weight_last, proj) overload for range safety")]]
double weighted_variance(Iterator first, Iterator last, WeightIterator weight_first, Projection proj)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::weighted_variance: empty range");
    }
 
    // Compute weighted mean
    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_variance: negative weight");
        }
 
        sum_weighted += value * weight;
        sum_weights += weight;
    }
 
    if (sum_weights == 0.0) {
        throw std::invalid_argument("statcpp::weighted_variance: sum of weights is zero");
    }
 
    double mean = sum_weighted / sum_weights;
 
    // Compute weighted variance
    double sum_squared_dev = 0.0;
    double sum_weights_squared = 0.0;
    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);
        double dev = value - mean;
 
        sum_squared_dev += weight * dev * dev;
        sum_weights_squared += weight * weight;
    }
 
    double correction = sum_weights - (sum_weights_squared / sum_weights);
 
    if (correction <= 0.0) {
        throw std::invalid_argument("statcpp::weighted_variance: insufficient effective sample size");
    }
 
    return sum_squared_dev / correction;
}
 
template <typename Iterator, typename WeightIterator>
double weighted_stddev(Iterator first, Iterator last, WeightIterator weight_first, WeightIterator weight_last)
{
    return std::sqrt(weighted_variance(first, last, weight_first, weight_last));
}
 
template <typename Iterator, typename WeightIterator>
[[deprecated("Use weighted_stddev(first, last, weight_first, weight_last) overload for range safety")]]
double weighted_stddev(Iterator first, Iterator last, WeightIterator weight_first)
{
    return std::sqrt(weighted_variance(first, last, weight_first));
}
 
template <typename Iterator, typename WeightIterator, typename Projection>
double weighted_stddev(Iterator first, Iterator last, WeightIterator weight_first, WeightIterator weight_last, Projection proj)
{
    return std::sqrt(weighted_variance(first, last, weight_first, weight_last, proj));
}
 
template <typename Iterator, typename WeightIterator, typename Projection>
[[deprecated("Use weighted_stddev(first, last, weight_first, weight_last, proj) overload for range safety")]]
double weighted_stddev(Iterator first, Iterator last, WeightIterator weight_first, Projection proj)
{
    return std::sqrt(weighted_variance(first, last, weight_first, proj));
}
 
// ============================================================================
// Geometric Standard Deviation
// ============================================================================
 
template <typename Iterator>
double geometric_stddev(Iterator first, Iterator last)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::geometric_stddev: empty range");
    }
 
    // Compute standard deviation of log-transformed data
    std::vector<double> log_values;
    log_values.reserve(n);
 
    for (auto it = first; it != last; ++it) {
        double value = static_cast<double>(*it);
        if (value <= 0.0) {
            throw std::invalid_argument("statcpp::geometric_stddev: all values must be positive");
        }
        log_values.push_back(std::log(value));
    }
 
    double log_stddev = statcpp::sample_stddev(log_values.begin(), log_values.end());
    return std::exp(log_stddev);
}
 
template <typename Iterator, typename Projection>
double geometric_stddev(Iterator first, Iterator last, Projection proj)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::geometric_stddev: empty range");
    }
 
    // Compute standard deviation of log-transformed data
    std::vector<double> log_values;
    log_values.reserve(n);
 
    for (auto it = first; it != last; ++it) {
        double value = static_cast<double>(std::invoke(proj, *it));
        if (value <= 0.0) {
            throw std::invalid_argument("statcpp::geometric_stddev: all values must be positive");
        }
        log_values.push_back(std::log(value));
    }
 
    double log_stddev = statcpp::sample_stddev(log_values.begin(), log_values.end());
    return std::exp(log_stddev);
}
 
} // namespace statcpp