estimation.hpp

Go to the documentation of this file.

 
#pragma once
 
#include "statcpp/basic_statistics.hpp"
#include "statcpp/dispersion_spread.hpp"
#include "statcpp/continuous_distributions.hpp"
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <stdexcept>
 
namespace statcpp {
 
// ============================================================================
// Standard Error
// ============================================================================
 
template <typename Iterator>
double standard_error(Iterator first, Iterator last)
{
    auto n = statcpp::count(first, last);
    if (n < 2) {
        throw std::invalid_argument("statcpp::standard_error: need at least 2 elements");
    }
    double s = statcpp::sample_stddev(first, last);
    return s / std::sqrt(static_cast<double>(n));
}
 
template <typename Iterator, typename Projection>
double standard_error(Iterator first, Iterator last, Projection proj)
{
    auto n = statcpp::count(first, last);
    if (n < 2) {
        throw std::invalid_argument("statcpp::standard_error: need at least 2 elements");
    }
    double s = statcpp::sample_stddev(first, last, proj);
    return s / std::sqrt(static_cast<double>(n));
}
 
template <typename Iterator>
double standard_error(Iterator first, Iterator last, double precomputed_stddev)
{
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::standard_error: empty range");
    }
    return precomputed_stddev / std::sqrt(static_cast<double>(n));
}
 
// ============================================================================
// Confidence Interval Result Structure
// ============================================================================
 
struct confidence_interval {
    double lower;              
    double upper;              
    double point_estimate;     
    double confidence_level;   
};
 
// ============================================================================
// Confidence Interval for Mean (using t-distribution)
// ============================================================================
 
template <typename Iterator>
confidence_interval ci_mean(Iterator first, Iterator last, double confidence = 0.95)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::ci_mean: confidence must be in (0, 1)");
    }
 
    auto n = statcpp::count(first, last);
    if (n < 2) {
        throw std::invalid_argument("statcpp::ci_mean: need at least 2 elements");
    }
 
    double mean_val = statcpp::mean(first, last);
    double se = statcpp::standard_error(first, last);
    double df = static_cast<double>(n - 1);
 
    double alpha = 1.0 - confidence;
    double t_crit = t_quantile(1.0 - alpha / 2.0, df);
 
    double margin = t_crit * se;
 
    return {mean_val - margin, mean_val + margin, mean_val, confidence};
}
 
template <typename Iterator, typename Projection>
confidence_interval ci_mean(Iterator first, Iterator last, double confidence, Projection proj)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::ci_mean: confidence must be in (0, 1)");
    }
 
    auto n = statcpp::count(first, last);
    if (n < 2) {
        throw std::invalid_argument("statcpp::ci_mean: need at least 2 elements");
    }
 
    double mean_val = statcpp::mean(first, last, proj);
    double se = statcpp::standard_error(first, last, proj);
    double df = static_cast<double>(n - 1);
 
    double alpha = 1.0 - confidence;
    double t_crit = t_quantile(1.0 - alpha / 2.0, df);
 
    double margin = t_crit * se;
 
    return {mean_val - margin, mean_val + margin, mean_val, confidence};
}
 
// ============================================================================
// Confidence Interval for Mean (using z-distribution, known variance)
// ============================================================================
 
template <typename Iterator>
confidence_interval ci_mean_z(Iterator first, Iterator last, double sigma, double confidence = 0.95)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::ci_mean_z: confidence must be in (0, 1)");
    }
    if (sigma <= 0.0) {
        throw std::invalid_argument("statcpp::ci_mean_z: sigma must be positive");
    }
 
    auto n = statcpp::count(first, last);
    if (n == 0) {
        throw std::invalid_argument("statcpp::ci_mean_z: empty range");
    }
 
    double mean_val = statcpp::mean(first, last);
    double se = sigma / std::sqrt(static_cast<double>(n));
 
    double alpha = 1.0 - confidence;
    double z_crit = norm_quantile(1.0 - alpha / 2.0);
 
    double margin = z_crit * se;
 
    return {mean_val - margin, mean_val + margin, mean_val, confidence};
}
 
// ============================================================================
// Confidence Interval for Proportion
// ============================================================================
 
inline confidence_interval ci_proportion(std::size_t successes, std::size_t trials, double confidence = 0.95)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::ci_proportion: confidence must be in (0, 1)");
    }
    if (trials == 0) {
        throw std::invalid_argument("statcpp::ci_proportion: trials must be positive");
    }
    if (successes > trials) {
        throw std::invalid_argument("statcpp::ci_proportion: successes cannot exceed trials");
    }
 
    double p_hat = static_cast<double>(successes) / static_cast<double>(trials);
    double se = std::sqrt(p_hat * (1.0 - p_hat) / static_cast<double>(trials));
 
    double alpha = 1.0 - confidence;
    double z_crit = norm_quantile(1.0 - alpha / 2.0);
 
    double margin = z_crit * se;
    double lower = std::max(0.0, p_hat - margin);
    double upper = std::min(1.0, p_hat + margin);
 
    return {lower, upper, p_hat, confidence};
}
 
inline confidence_interval ci_proportion_wilson(std::size_t successes, std::size_t trials, double confidence = 0.95)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::ci_proportion_wilson: confidence must be in (0, 1)");
    }
    if (trials == 0) {
        throw std::invalid_argument("statcpp::ci_proportion_wilson: trials must be positive");
    }
    if (successes > trials) {
        throw std::invalid_argument("statcpp::ci_proportion_wilson: successes cannot exceed trials");
    }
 
    double n = static_cast<double>(trials);
    double p_hat = static_cast<double>(successes) / n;
 
    double alpha = 1.0 - confidence;
    double z = norm_quantile(1.0 - alpha / 2.0);
    double z2 = z * z;
 
    double denom = 1.0 + z2 / n;
    double center = (p_hat + z2 / (2.0 * n)) / denom;
    double margin = z * std::sqrt((p_hat * (1.0 - p_hat) + z2 / (4.0 * n)) / n) / denom;
 
    return {center - margin, center + margin, p_hat, confidence};
}
 
// ============================================================================
// Confidence Interval for Variance (Chi-square based)
// ============================================================================
 
template <typename Iterator>
confidence_interval ci_variance(Iterator first, Iterator last, double confidence = 0.95)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::ci_variance: confidence must be in (0, 1)");
    }
 
    auto n = statcpp::count(first, last);
    if (n < 2) {
        throw std::invalid_argument("statcpp::ci_variance: need at least 2 elements");
    }
 
    double var = statcpp::sample_variance(first, last);
    double df = static_cast<double>(n - 1);
 
    double alpha = 1.0 - confidence;
    double chi2_lower = chisq_quantile(alpha / 2.0, df);
    double chi2_upper = chisq_quantile(1.0 - alpha / 2.0, df);
 
    double lower = df * var / chi2_upper;
    double upper = df * var / chi2_lower;
 
    return {lower, upper, var, confidence};
}
 
// ============================================================================
// Confidence Interval for Difference of Means
// ============================================================================
 
template <typename Iterator1, typename Iterator2>
confidence_interval ci_mean_diff(Iterator1 first1, Iterator1 last1,
                                  Iterator2 first2, Iterator2 last2,
                                  double confidence = 0.95)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::ci_mean_diff: confidence must be in (0, 1)");
    }
 
    auto n1 = statcpp::count(first1, last1);
    auto n2 = statcpp::count(first2, last2);
 
    if (n1 < 2 || n2 < 2) {
        throw std::invalid_argument("statcpp::ci_mean_diff: need at least 2 elements in each sample");
    }
 
    double mean1 = statcpp::mean(first1, last1);
    double mean2 = statcpp::mean(first2, last2);
    double var1 = statcpp::sample_variance(first1, last1);
    double var2 = statcpp::sample_variance(first2, last2);
 
    double diff = mean1 - mean2;
 
    // Pooled variance
    double df = static_cast<double>(n1 + n2 - 2);
    double sp2 = ((n1 - 1) * var1 + (n2 - 1) * var2) / df;
    double se = std::sqrt(sp2 * (1.0 / n1 + 1.0 / n2));
 
    double alpha = 1.0 - confidence;
    double t_crit = t_quantile(1.0 - alpha / 2.0, df);
 
    double margin = t_crit * se;
 
    return {diff - margin, diff + margin, diff, confidence};
}
 
template <typename Iterator1, typename Iterator2>
confidence_interval ci_mean_diff_welch(Iterator1 first1, Iterator1 last1,
                                        Iterator2 first2, Iterator2 last2,
                                        double confidence = 0.95)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::ci_mean_diff_welch: confidence must be in (0, 1)");
    }
 
    auto n1 = statcpp::count(first1, last1);
    auto n2 = statcpp::count(first2, last2);
 
    if (n1 < 2 || n2 < 2) {
        throw std::invalid_argument("statcpp::ci_mean_diff_welch: need at least 2 elements in each sample");
    }
 
    double mean1 = statcpp::mean(first1, last1);
    double mean2 = statcpp::mean(first2, last2);
    double var1 = statcpp::sample_variance(first1, last1);
    double var2 = statcpp::sample_variance(first2, last2);
 
    double diff = mean1 - mean2;
 
    double se1 = var1 / n1;
    double se2 = var2 / n2;
    double se = std::sqrt(se1 + se2);
 
    // Welch-Satterthwaite approximation for degrees of freedom
    double num = (se1 + se2) * (se1 + se2);
    double denom = (se1 * se1) / (n1 - 1) + (se2 * se2) / (n2 - 1);
 
    // Protection against zero denominator (when both variances are zero)
    if (denom == 0.0) {
        throw std::invalid_argument("statcpp::ci_mean_diff_welch: cannot compute degrees of freedom with zero variances");
    }
 
    double df = num / denom;
 
    double alpha = 1.0 - confidence;
    double t_crit = t_quantile(1.0 - alpha / 2.0, df);
 
    double margin = t_crit * se;
 
    return {diff - margin, diff + margin, diff, confidence};
}
 
// ============================================================================
// Margin of Error
// ============================================================================
 
template <typename Iterator>
double margin_of_error_mean(Iterator first, Iterator last, double confidence = 0.95)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::margin_of_error_mean: confidence must be in (0, 1)");
    }
 
    auto n = statcpp::count(first, last);
    if (n < 2) {
        throw std::invalid_argument("statcpp::margin_of_error_mean: need at least 2 elements");
    }
 
    double se = statcpp::standard_error(first, last);
    double df = static_cast<double>(n - 1);
    double alpha = 1.0 - confidence;
    double t_crit = t_quantile(1.0 - alpha / 2.0, df);
 
    return t_crit * se;
}
 
template <typename Iterator, typename Projection>
double margin_of_error_mean(Iterator first, Iterator last, double confidence, Projection proj)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::margin_of_error_mean: confidence must be in (0, 1)");
    }
 
    auto n = statcpp::count(first, last);
    if (n < 2) {
        throw std::invalid_argument("statcpp::margin_of_error_mean: need at least 2 elements");
    }
 
    double se = statcpp::standard_error(first, last, proj);
    double df = static_cast<double>(n - 1);
    double alpha = 1.0 - confidence;
    double t_crit = t_quantile(1.0 - alpha / 2.0, df);
 
    return t_crit * se;
}
 
inline double margin_of_error_proportion(std::size_t successes, std::size_t n, double confidence = 0.95)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::margin_of_error_proportion: confidence must be in (0, 1)");
    }
    if (n == 0) {
        throw std::invalid_argument("statcpp::margin_of_error_proportion: n must be positive");
    }
    if (successes > n) {
        throw std::invalid_argument("statcpp::margin_of_error_proportion: successes cannot exceed n");
    }
 
    double p = static_cast<double>(successes) / static_cast<double>(n);
    double se = std::sqrt(p * (1.0 - p) / static_cast<double>(n));
 
    double alpha = 1.0 - confidence;
    double z_crit = normal_quantile(1.0 - alpha / 2.0, 0.0, 1.0);
 
    return z_crit * se;
}
 
inline double margin_of_error_proportion_worst_case(std::size_t n, double confidence = 0.95)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::margin_of_error_proportion_worst_case: confidence must be in (0, 1)");
    }
    if (n == 0) {
        throw std::invalid_argument("statcpp::margin_of_error_proportion_worst_case: n must be positive");
    }
 
    double alpha = 1.0 - confidence;
    double z_crit = normal_quantile(1.0 - alpha / 2.0, 0.0, 1.0);
 
    return z_crit * 0.5 / std::sqrt(static_cast<double>(n));
}
 
// ============================================================================
// Sample Size Calculation (Margin-of-Error based)
// ============================================================================
 
inline std::size_t sample_size_for_moe_proportion(double margin_of_error,
                                                   double confidence_level = 0.95,
                                                   double p_estimate = 0.5)
{
    if (margin_of_error <= 0.0 || margin_of_error >= 1.0) {
        throw std::invalid_argument("statcpp::sample_size_for_moe_proportion: margin_of_error must be in (0, 1)");
    }
    if (confidence_level <= 0.0 || confidence_level >= 1.0) {
        throw std::invalid_argument("statcpp::sample_size_for_moe_proportion: confidence_level must be in (0, 1)");
    }
    if (p_estimate <= 0.0 || p_estimate >= 1.0) {
        throw std::invalid_argument("statcpp::sample_size_for_moe_proportion: p_estimate must be in (0, 1)");
    }
 
    double alpha = 1.0 - confidence_level;
    double z = normal_quantile(1.0 - alpha / 2.0, 0.0, 1.0);
 
    // n = (z / MoE)^2 * p(1-p)
    double n_exact = std::pow(z / margin_of_error, 2.0) * p_estimate * (1.0 - p_estimate);
 
    // Round up to integer
    return static_cast<std::size_t>(std::ceil(n_exact));
}
 
inline std::size_t sample_size_for_moe_mean(double margin_of_error,
                                             double sigma,
                                             double confidence_level = 0.95)
{
    if (margin_of_error <= 0.0) {
        throw std::invalid_argument("statcpp::sample_size_for_moe_mean: margin_of_error must be positive");
    }
    if (sigma <= 0.0) {
        throw std::invalid_argument("statcpp::sample_size_for_moe_mean: sigma must be positive");
    }
    if (confidence_level <= 0.0 || confidence_level >= 1.0) {
        throw std::invalid_argument("statcpp::sample_size_for_moe_mean: confidence_level must be in (0, 1)");
    }
 
    double alpha = 1.0 - confidence_level;
    double z = normal_quantile(1.0 - alpha / 2.0, 0.0, 1.0);
 
    // n = (z * sigma / MoE)^2
    double n_exact = std::pow(z * sigma / margin_of_error, 2.0);
 
    // Round up to integer
    return static_cast<std::size_t>(std::ceil(n_exact));
}
 
// ============================================================================
// Two-Sample Mean Difference Confidence Interval
// ============================================================================
 
template <typename Iterator1, typename Iterator2>
confidence_interval ci_mean_diff_pooled(Iterator1 first1, Iterator1 last1,
                                        Iterator2 first2, Iterator2 last2,
                                        double confidence = 0.95)
{
    // ci_mean_diff と同一のロジック(等分散を仮定した pooled t 検定)
    return ci_mean_diff(first1, last1, first2, last2, confidence);
}
 
// ============================================================================
// Two-Sample Proportion Difference Confidence Interval
// ============================================================================
 
inline confidence_interval ci_proportion_diff(std::size_t successes1, std::size_t n1,
                                               std::size_t successes2, std::size_t n2,
                                               double confidence = 0.95)
{
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("statcpp::ci_proportion_diff: confidence must be in (0, 1)");
    }
    if (n1 == 0 || n2 == 0) {
        throw std::invalid_argument("statcpp::ci_proportion_diff: sample sizes must be positive");
    }
    if (successes1 > n1 || successes2 > n2) {
        throw std::invalid_argument("statcpp::ci_proportion_diff: successes cannot exceed sample size");
    }
 
    double p1 = static_cast<double>(successes1) / static_cast<double>(n1);
    double p2 = static_cast<double>(successes2) / static_cast<double>(n2);
 
    double diff = p1 - p2;
 
    // Standard error for difference of proportions
    double se1 = p1 * (1.0 - p1) / static_cast<double>(n1);
    double se2 = p2 * (1.0 - p2) / static_cast<double>(n2);
    double se = std::sqrt(se1 + se2);
 
    double alpha = 1.0 - confidence;
    double z_crit = normal_quantile(1.0 - alpha / 2.0, 0.0, 1.0);
 
    double margin = z_crit * se;
 
    return {diff - margin, diff + margin, diff, confidence};
}
 
} // namespace statcpp