robust.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <iterator>
#include <limits>
#include <stdexcept>
#include <vector>
 
#include "statcpp/basic_statistics.hpp"
#include "statcpp/dispersion_spread.hpp"
#include "statcpp/order_statistics.hpp"
 
namespace statcpp {
 
// ============================================================================
// Median Absolute Deviation (MAD)
// ============================================================================
 
template <typename Iterator>
double mad(Iterator first, Iterator last)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::mad: empty range");
    }
 
    // Copy and sort data (for median calculation)
    std::vector<double> sorted_data;
    sorted_data.reserve(n);
    for (auto it = first; it != last; ++it) {
        sorted_data.push_back(static_cast<double>(*it));
    }
    std::sort(sorted_data.begin(), sorted_data.end());
 
    double med = statcpp::median(sorted_data.begin(), sorted_data.end());
 
    // Calculate absolute deviations from median
    std::vector<double> abs_deviations;
    abs_deviations.reserve(n);
    for (double val : sorted_data) {
        abs_deviations.push_back(std::abs(val - med));
    }
 
    // Median of absolute deviations
    std::sort(abs_deviations.begin(), abs_deviations.end());
    return statcpp::median(abs_deviations.begin(), abs_deviations.end());
}
 
template <typename Iterator>
double mad_scaled(Iterator first, Iterator last)
{
    return 1.4826 * mad(first, last);
}
 
// ============================================================================
// Outlier Detection (IQR Method / Tukey's Fences)
// ============================================================================
 
struct outlier_detection_result {
    std::vector<double> outliers;                  
    std::vector<std::size_t> outlier_indices;      
    double lower_fence;                            
    double upper_fence;                            
    double q1;                                     
    double q3;                                     
    double iqr_value;                              
};
 
template <typename Iterator>
outlier_detection_result detect_outliers_iqr(Iterator first, Iterator last, double k = 1.5)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::detect_outliers_iqr: empty range");
    }
 
    // Copy and sort data
    std::vector<double> sorted_data;
    sorted_data.reserve(n);
    for (auto it = first; it != last; ++it) {
        sorted_data.push_back(static_cast<double>(*it));
    }
    std::sort(sorted_data.begin(), sorted_data.end());
 
    auto q = statcpp::quartiles(sorted_data.begin(), sorted_data.end());
    double iqr_val = q.q3 - q.q1;
 
    double lower_fence = q.q1 - k * iqr_val;
    double upper_fence = q.q3 + k * iqr_val;
 
    outlier_detection_result result;
    result.lower_fence = lower_fence;
    result.upper_fence = upper_fence;
    result.q1 = q.q1;
    result.q3 = q.q3;
    result.iqr_value = iqr_val;
 
    // Detect outliers
    std::size_t idx = 0;
    for (auto it = first; it != last; ++it, ++idx) {
        double val = static_cast<double>(*it);
        if (val < lower_fence || val > upper_fence) {
            result.outliers.push_back(val);
            result.outlier_indices.push_back(idx);
        }
    }
 
    return result;
}
 
// ============================================================================
// Z-score Outlier Detection
// ============================================================================
 
template <typename Iterator>
outlier_detection_result detect_outliers_zscore(Iterator first, Iterator last, double threshold = 3.0)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n < 2) {
        throw std::invalid_argument("statcpp::detect_outliers_zscore: need at least 2 elements");
    }
 
    double m = statcpp::mean(first, last);
    double s = statcpp::sample_stddev(first, last, m);
 
    if (s == 0.0) {
        throw std::invalid_argument("statcpp::detect_outliers_zscore: zero standard deviation");
    }
 
    outlier_detection_result result;
    result.lower_fence = m - threshold * s;
    result.upper_fence = m + threshold * s;
    result.q1 = 0.0;  // Z-score method does not use quartiles
    result.q3 = 0.0;
    result.iqr_value = 0.0;
 
    std::size_t idx = 0;
    for (auto it = first; it != last; ++it, ++idx) {
        double val = static_cast<double>(*it);
        double z = (val - m) / s;
        if (std::abs(z) > threshold) {
            result.outliers.push_back(val);
            result.outlier_indices.push_back(idx);
        }
    }
 
    return result;
}
 
template <typename Iterator>
outlier_detection_result detect_outliers_modified_zscore(Iterator first, Iterator last, double threshold = 3.5)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::detect_outliers_modified_zscore: empty range");
    }
 
    // Copy and sort data
    std::vector<double> sorted_data;
    sorted_data.reserve(n);
    for (auto it = first; it != last; ++it) {
        sorted_data.push_back(static_cast<double>(*it));
    }
    std::sort(sorted_data.begin(), sorted_data.end());
 
    double med = statcpp::median(sorted_data.begin(), sorted_data.end());
    double mad_val = mad(first, last);
 
    if (mad_val == 0.0) {
        throw std::invalid_argument("statcpp::detect_outliers_modified_zscore: zero MAD");
    }
 
    // Modified Z-score = 0.6745 * (x - median) / MAD
    double scale = 0.6745;
 
    outlier_detection_result result;
    result.lower_fence = med - threshold * mad_val / scale;
    result.upper_fence = med + threshold * mad_val / scale;
    result.q1 = 0.0;
    result.q3 = 0.0;
    result.iqr_value = 0.0;
 
    std::size_t idx = 0;
    for (auto it = first; it != last; ++it, ++idx) {
        double val = static_cast<double>(*it);
        double modified_z = scale * (val - med) / mad_val;
        if (std::abs(modified_z) > threshold) {
            result.outliers.push_back(val);
            result.outlier_indices.push_back(idx);
        }
    }
 
    return result;
}
 
// ============================================================================
// Winsorization
// ============================================================================
 
template <typename Iterator>
std::vector<double> winsorize(Iterator first, Iterator last, double limits = 0.05)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::winsorize: empty range");
    }
    if (limits < 0.0 || limits >= 0.5) {
        throw std::invalid_argument("statcpp::winsorize: limits must be in [0, 0.5)");
    }
 
    // Copy and sort data
    std::vector<double> sorted_data;
    sorted_data.reserve(n);
    for (auto it = first; it != last; ++it) {
        sorted_data.push_back(static_cast<double>(*it));
    }
    std::sort(sorted_data.begin(), sorted_data.end());
 
    // Calculate thresholds
    double lower_percentile = limits;
    double upper_percentile = 1.0 - limits;
    double lower_value = statcpp::percentile(sorted_data.begin(), sorted_data.end(), lower_percentile);
    double upper_value = statcpp::percentile(sorted_data.begin(), sorted_data.end(), upper_percentile);
 
    // Winsorize
    std::vector<double> result;
    result.reserve(n);
    for (auto it = first; it != last; ++it) {
        double val = static_cast<double>(*it);
        if (val < lower_value) {
            result.push_back(lower_value);
        } else if (val > upper_value) {
            result.push_back(upper_value);
        } else {
            result.push_back(val);
        }
    }
 
    return result;
}
 
// ============================================================================
// Cook's Distance (for Linear Regression Diagnostics)
// ============================================================================
 
inline std::vector<double> cooks_distance(
    const std::vector<double>& residuals,
    const std::vector<double>& hat_values,
    double mse,
    std::size_t p)
{
    if (residuals.size() != hat_values.size()) {
        throw std::invalid_argument("statcpp::cooks_distance: residuals and hat_values must have same length");
    }
    if (residuals.empty()) {
        throw std::invalid_argument("statcpp::cooks_distance: empty data");
    }
    if (mse <= 0) {
        throw std::invalid_argument("statcpp::cooks_distance: mse must be positive");
    }
    if (p == 0) {
        throw std::invalid_argument("statcpp::cooks_distance: p must be positive");
    }
 
    std::vector<double> result;
    result.reserve(residuals.size());
 
    for (std::size_t i = 0; i < residuals.size(); ++i) {
        double h = hat_values[i];
        if (h >= 1.0) {
            result.push_back(std::numeric_limits<double>::infinity());
            continue;
        }
 
        double e = residuals[i];
        double d = (e * e / (static_cast<double>(p) * mse)) *
                   (h / ((1.0 - h) * (1.0 - h)));
        result.push_back(d);
    }
 
    return result;
}
 
// ============================================================================
// DFFITS (Difference in Fits)
// ============================================================================
 
inline std::vector<double> dffits(
    const std::vector<double>& residuals,
    const std::vector<double>& hat_values,
    double mse)
{
    if (residuals.size() != hat_values.size()) {
        throw std::invalid_argument("statcpp::dffits: residuals and hat_values must have same length");
    }
    if (residuals.empty()) {
        throw std::invalid_argument("statcpp::dffits: empty data");
    }
 
    std::vector<double> result;
    result.reserve(residuals.size());
 
    for (std::size_t i = 0; i < residuals.size(); ++i) {
        double h = hat_values[i];
        if (h >= 1.0) {
            result.push_back(std::numeric_limits<double>::infinity());
            continue;
        }
 
        double e = residuals[i];
        // Studentized residual
        double se = std::sqrt(mse * (1.0 - h));
        if (se == 0.0) {
            result.push_back(std::numeric_limits<double>::infinity());
            continue;
        }
        double t = e / se;
 
        // DFFITS
        double dffits_val = t * std::sqrt(h / (1.0 - h));
        result.push_back(dffits_val);
    }
 
    return result;
}
 
// ============================================================================
// Robust Location Estimators
// ============================================================================
 
template <typename Iterator>
double hodges_lehmann(Iterator first, Iterator last)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n == 0) {
        throw std::invalid_argument("statcpp::hodges_lehmann: empty range");
    }
 
    std::vector<double> data;
    data.reserve(n);
    for (auto it = first; it != last; ++it) {
        data.push_back(static_cast<double>(*it));
    }
 
    // Calculate all pairwise averages
    std::vector<double> pairwise_means;
    pairwise_means.reserve(n * (n + 1) / 2);
 
    for (std::size_t i = 0; i < n; ++i) {
        for (std::size_t j = i; j < n; ++j) {
            pairwise_means.push_back((data[i] + data[j]) / 2.0);
        }
    }
 
    std::sort(pairwise_means.begin(), pairwise_means.end());
    return statcpp::median(pairwise_means.begin(), pairwise_means.end());
}
 
template <typename Iterator>
double biweight_midvariance(Iterator first, Iterator last, double c = 9.0)
{
    auto n = static_cast<std::size_t>(std::distance(first, last));
    if (n < 2) {
        throw std::invalid_argument("statcpp::biweight_midvariance: need at least 2 elements");
    }
 
    // Copy and sort data
    std::vector<double> data;
    data.reserve(n);
    for (auto it = first; it != last; ++it) {
        data.push_back(static_cast<double>(*it));
    }
    std::sort(data.begin(), data.end());
 
    double med = statcpp::median(data.begin(), data.end());
    double mad_val = mad(first, last);
 
    if (mad_val == 0.0) {
        return 0.0;
    }
 
    double num = 0.0;
    double den = 0.0;
    std::size_t count = 0;
 
    for (double x : data) {
        double u = (x - med) / (c * mad_val);
        if (std::abs(u) < 1.0) {
            double u2 = u * u;
            double w = (1.0 - u2);
            double w2 = w * w;
            double w4 = w2 * w2;
            num += (x - med) * (x - med) * w4;
            den += w * (1.0 - 5.0 * u2);
            count++;
        }
    }
 
    if (den == 0.0 || count < 2) {
        return 0.0;
    }
 
    return static_cast<double>(n) * num / (den * den);
}
 
} // namespace statcpp