continuous_distributions.hpp

Go to the documentation of this file.

 
#pragma once
 
#include "statcpp/special_functions.hpp"
#include "statcpp/random_engine.hpp"
 
#include <algorithm>
#include <cmath>
#include <limits>
#include <random>
#include <stdexcept>
 
namespace statcpp {
 
// ============================================================================
// Uniform Distribution
// ============================================================================
 
inline double uniform_pdf(double x, double a = 0.0, double b = 1.0)
{
    if (a >= b) {
        throw std::invalid_argument("statcpp::uniform_pdf: a must be less than b");
    }
    if (x < a || x > b) {
        return 0.0;
    }
    return 1.0 / (b - a);
}
 
inline double uniform_cdf(double x, double a = 0.0, double b = 1.0)
{
    if (a >= b) {
        throw std::invalid_argument("statcpp::uniform_cdf: a must be less than b");
    }
    if (x < a) return 0.0;
    if (x > b) return 1.0;
    return (x - a) / (b - a);
}
 
inline double uniform_quantile(double p, double a = 0.0, double b = 1.0)
{
    if (a >= b) {
        throw std::invalid_argument("statcpp::uniform_quantile: a must be less than b");
    }
    if (p < 0.0 || p > 1.0) {
        throw std::invalid_argument("statcpp::uniform_quantile: p must be in [0, 1]");
    }
    return a + p * (b - a);
}
 
template <typename Engine = default_random_engine>
double uniform_rand(double a, double b, Engine& engine)
{
    if (a >= b) {
        throw std::invalid_argument("statcpp::uniform_rand: a must be less than b");
    }
    std::uniform_real_distribution<double> dist(a, b);
    return dist(engine);
}
 
inline double uniform_rand(double a = 0.0, double b = 1.0)
{
    return uniform_rand(a, b, get_random_engine());
}
 
// ============================================================================
// Normal Distribution
// ============================================================================
 
inline double normal_pdf(double x, double mu = 0.0, double sigma = 1.0)
{
    if (sigma <= 0.0) {
        throw std::invalid_argument("statcpp::normal_pdf: sigma must be positive");
    }
    double z = (x - mu) / sigma;
    return std::exp(-0.5 * z * z) / (sigma * sqrt_2_pi);
}
 
inline double normal_cdf(double x, double mu = 0.0, double sigma = 1.0)
{
    if (sigma <= 0.0) {
        throw std::invalid_argument("statcpp::normal_cdf: sigma must be positive");
    }
    return norm_cdf((x - mu) / sigma);
}
 
inline double normal_quantile(double p, double mu = 0.0, double sigma = 1.0)
{
    if (sigma <= 0.0) {
        throw std::invalid_argument("statcpp::normal_quantile: sigma must be positive");
    }
    if (p <= 0.0 || p >= 1.0) {
        if (p == 0.0) return -std::numeric_limits<double>::infinity();
        if (p == 1.0) return std::numeric_limits<double>::infinity();
        throw std::invalid_argument("statcpp::normal_quantile: p must be in (0, 1)");
    }
    return mu + sigma * norm_quantile(p);
}
 
template <typename Engine = default_random_engine>
double normal_rand(double mu, double sigma, Engine& engine)
{
    if (sigma <= 0.0) {
        throw std::invalid_argument("statcpp::normal_rand: sigma must be positive");
    }
    std::normal_distribution<double> dist(mu, sigma);
    return dist(engine);
}
 
inline double normal_rand(double mu = 0.0, double sigma = 1.0)
{
    return normal_rand(mu, sigma, get_random_engine());
}
 
// ============================================================================
// Exponential Distribution
// ============================================================================
 
inline double exponential_pdf(double x, double lambda = 1.0)
{
    if (lambda <= 0.0) {
        throw std::invalid_argument("statcpp::exponential_pdf: lambda must be positive");
    }
    if (x < 0.0) return 0.0;
    return lambda * std::exp(-lambda * x);
}
 
inline double exponential_cdf(double x, double lambda = 1.0)
{
    if (lambda <= 0.0) {
        throw std::invalid_argument("statcpp::exponential_cdf: lambda must be positive");
    }
    if (x < 0.0) return 0.0;
    return -std::expm1(-lambda * x);
}
 
inline double exponential_quantile(double p, double lambda = 1.0)
{
    if (lambda <= 0.0) {
        throw std::invalid_argument("statcpp::exponential_quantile: lambda must be positive");
    }
    if (p < 0.0 || p >= 1.0) {
        if (p == 1.0) return std::numeric_limits<double>::infinity();
        throw std::invalid_argument("statcpp::exponential_quantile: p must be in [0, 1)");
    }
    return -std::log(1.0 - p) / lambda;
}
 
template <typename Engine = default_random_engine>
double exponential_rand(double lambda, Engine& engine)
{
    if (lambda <= 0.0) {
        throw std::invalid_argument("statcpp::exponential_rand: lambda must be positive");
    }
    std::exponential_distribution<double> dist(lambda);
    return dist(engine);
}
 
inline double exponential_rand(double lambda = 1.0)
{
    return exponential_rand(lambda, get_random_engine());
}
 
// ============================================================================
// Gamma Distribution
// ============================================================================
 
inline double gamma_pdf(double x, double shape, double rate = 1.0)
{
    if (shape <= 0.0) {
        throw std::invalid_argument("statcpp::gamma_pdf: shape must be positive");
    }
    if (rate <= 0.0) {
        throw std::invalid_argument("statcpp::gamma_pdf: rate must be positive");
    }
    if (x < 0.0) return 0.0;
    if (x == 0.0) {
        if (shape < 1.0) return std::numeric_limits<double>::infinity();
        if (shape == 1.0) return rate;  // Exponential distribution PDF at x=0
        return 0.0;  // shape > 1
    }
 
    return std::exp(shape * std::log(rate) + (shape - 1.0) * std::log(x) - rate * x - lgamma(shape));
}
 
inline double gamma_cdf(double x, double shape, double rate = 1.0)
{
    if (shape <= 0.0) {
        throw std::invalid_argument("statcpp::gamma_cdf: shape must be positive");
    }
    if (rate <= 0.0) {
        throw std::invalid_argument("statcpp::gamma_cdf: rate must be positive");
    }
    if (x <= 0.0) return 0.0;
 
    return gammainc_lower(shape, rate * x);
}
 
inline double gamma_quantile(double p, double shape, double rate = 1.0)
{
    if (shape <= 0.0) {
        throw std::invalid_argument("statcpp::gamma_quantile: shape must be positive");
    }
    if (rate <= 0.0) {
        throw std::invalid_argument("statcpp::gamma_quantile: rate must be positive");
    }
    if (p < 0.0 || p > 1.0) {
        throw std::invalid_argument("statcpp::gamma_quantile: p must be in [0, 1]");
    }
    if (p == 0.0) return 0.0;
    if (p == 1.0) return std::numeric_limits<double>::infinity();
 
    return gammainc_lower_inv(shape, p) / rate;
}
 
template <typename Engine = default_random_engine>
double gamma_rand(double shape, double rate, Engine& engine)
{
    if (shape <= 0.0) {
        throw std::invalid_argument("statcpp::gamma_rand: shape must be positive");
    }
    if (rate <= 0.0) {
        throw std::invalid_argument("statcpp::gamma_rand: rate must be positive");
    }
    // std::gamma_distribution uses shape and scale (1/rate)
    std::gamma_distribution<double> dist(shape, 1.0 / rate);
    return dist(engine);
}
 
inline double gamma_rand(double shape, double rate = 1.0)
{
    return gamma_rand(shape, rate, get_random_engine());
}
 
// ============================================================================
// Beta Distribution
// ============================================================================
 
inline double beta_pdf(double x, double alpha, double beta_param)
{
    if (alpha <= 0.0) {
        throw std::invalid_argument("statcpp::beta_pdf: alpha must be positive");
    }
    if (beta_param <= 0.0) {
        throw std::invalid_argument("statcpp::beta_pdf: beta must be positive");
    }
    if (x < 0.0 || x > 1.0) return 0.0;
    if (x == 0.0) {
        if (alpha < 1.0) return std::numeric_limits<double>::infinity();
        if (alpha == 1.0) return beta_param;  // Beta(1, beta) PDF at x=0
        return 0.0;  // alpha > 1
    }
    if (x == 1.0) {
        if (beta_param < 1.0) return std::numeric_limits<double>::infinity();
        if (beta_param == 1.0) return alpha;  // Beta(alpha, 1) PDF at x=1
        return 0.0;  // beta > 1
    }
 
    return std::exp((alpha - 1.0) * std::log(x) + (beta_param - 1.0) * std::log(1.0 - x) - lbeta(alpha, beta_param));
}
 
inline double beta_cdf(double x, double alpha, double beta_param)
{
    if (alpha <= 0.0) {
        throw std::invalid_argument("statcpp::beta_cdf: alpha must be positive");
    }
    if (beta_param <= 0.0) {
        throw std::invalid_argument("statcpp::beta_cdf: beta must be positive");
    }
    if (x <= 0.0) return 0.0;
    if (x >= 1.0) return 1.0;
 
    return betainc(alpha, beta_param, x);
}
 
inline double beta_quantile(double p, double alpha, double beta_param)
{
    if (alpha <= 0.0) {
        throw std::invalid_argument("statcpp::beta_quantile: alpha must be positive");
    }
    if (beta_param <= 0.0) {
        throw std::invalid_argument("statcpp::beta_quantile: beta must be positive");
    }
    if (p < 0.0 || p > 1.0) {
        throw std::invalid_argument("statcpp::beta_quantile: p must be in [0, 1]");
    }
    if (p == 0.0) return 0.0;
    if (p == 1.0) return 1.0;
 
    return betaincinv(alpha, beta_param, p);
}
 
template <typename Engine = default_random_engine>
double beta_rand(double alpha, double beta_param, Engine& engine)
{
    if (alpha <= 0.0) {
        throw std::invalid_argument("statcpp::beta_rand: alpha must be positive");
    }
    if (beta_param <= 0.0) {
        throw std::invalid_argument("statcpp::beta_rand: beta must be positive");
    }
 
    std::gamma_distribution<double> dist_a(alpha, 1.0);
    std::gamma_distribution<double> dist_b(beta_param, 1.0);
 
    double x = dist_a(engine);
    double y = dist_b(engine);
    // Extremely small shape parameters can underflow both gamma variates to 0,
    // yielding 0/0 = NaN. Re-draw until the sum is positive.
    while (x + y <= 0.0) {
        x = dist_a(engine);
        y = dist_b(engine);
    }
    return x / (x + y);
}
 
inline double beta_rand(double alpha, double beta_param)
{
    return beta_rand(alpha, beta_param, get_random_engine());
}
 
// ============================================================================
// Chi-Square Distribution
// ============================================================================
 
inline double chisq_pdf(double x, double df)
{
    if (df <= 0.0) {
        throw std::invalid_argument("statcpp::chisq_pdf: df must be positive");
    }
    return gamma_pdf(x, df / 2.0, 0.5);
}
 
inline double chisq_cdf(double x, double df)
{
    if (df <= 0.0) {
        throw std::invalid_argument("statcpp::chisq_cdf: df must be positive");
    }
    return gamma_cdf(x, df / 2.0, 0.5);
}
 
inline double chisq_quantile(double p, double df)
{
    if (df <= 0.0) {
        throw std::invalid_argument("statcpp::chisq_quantile: df must be positive");
    }
    return gamma_quantile(p, df / 2.0, 0.5);
}
 
template <typename Engine = default_random_engine>
double chisq_rand(double df, Engine& engine)
{
    if (df <= 0.0) {
        throw std::invalid_argument("statcpp::chisq_rand: df must be positive");
    }
    return gamma_rand(df / 2.0, 0.5, engine);
}
 
inline double chisq_rand(double df)
{
    return chisq_rand(df, get_random_engine());
}
 
// ============================================================================
// Student's t-Distribution
// ============================================================================
 
inline double t_pdf(double x, double df)
{
    if (df <= 0.0) {
        throw std::invalid_argument("statcpp::t_pdf: df must be positive");
    }
 
    double coef = std::exp(lgamma((df + 1.0) / 2.0) - lgamma(df / 2.0)) / std::sqrt(df * pi);
    return coef * std::pow(1.0 + x * x / df, -(df + 1.0) / 2.0);
}
 
inline double t_cdf(double x, double df)
{
    if (df <= 0.0) {
        throw std::invalid_argument("statcpp::t_cdf: df must be positive");
    }
 
    double t2 = x * x;
    double p = betainc(df / 2.0, 0.5, df / (df + t2));
 
    if (x >= 0.0) {
        return 1.0 - 0.5 * p;
    } else {
        return 0.5 * p;
    }
}
 
inline double t_quantile(double p, double df)
{
    if (df <= 0.0) {
        throw std::invalid_argument("statcpp::t_quantile: df must be positive");
    }
    if (p <= 0.0 || p >= 1.0) {
        if (p == 0.0) return -std::numeric_limits<double>::infinity();
        if (p == 1.0) return std::numeric_limits<double>::infinity();
        throw std::invalid_argument("statcpp::t_quantile: p must be in (0, 1)");
    }
 
    // Use normal quantile as initial guess for large df
    double x = norm_quantile(p);
 
    // For small df (but df > 2), adjust initial guess using variance scaling
    if (df > 2.0 && df < 4.0) {
        x *= std::sqrt(df / (df - 2.0));
    }
 
    const double eps = 1e-10;
    const int max_iter = 50;
 
    for (int i = 0; i < max_iter; ++i) {
        if (!std::isfinite(x)) {
            x = norm_quantile(p);
        }
        double f = t_cdf(x, df) - p;
        if (std::abs(f) < eps) {
            return x;
        }
        double fprime = t_pdf(x, df);
        if (fprime == 0.0) break;
 
        double x_new = x - f / fprime;
 
        if (std::abs(x_new - x) < eps * (1.0 + std::abs(x))) {
            return x_new;
        }
 
        x = x_new;
    }
 
    if (!std::isfinite(x)) {
        x = norm_quantile(p);
    }
    return x;
}
 
template <typename Engine = default_random_engine>
double t_rand(double df, Engine& engine)
{
    if (df <= 0.0) {
        throw std::invalid_argument("statcpp::t_rand: df must be positive");
    }
    std::student_t_distribution<double> dist(df);
    return dist(engine);
}
 
inline double t_rand(double df)
{
    return t_rand(df, get_random_engine());
}
 
// ============================================================================
// F-Distribution
// ============================================================================
 
inline double f_pdf(double x, double df1, double df2)
{
    if (df1 <= 0.0) {
        throw std::invalid_argument("statcpp::f_pdf: df1 must be positive");
    }
    if (df2 <= 0.0) {
        throw std::invalid_argument("statcpp::f_pdf: df2 must be positive");
    }
    if (x <= 0.0) return 0.0;
 
    double log_pdf = (df1 / 2.0) * std::log(df1) + (df2 / 2.0) * std::log(df2)
                   + (df1 / 2.0 - 1.0) * std::log(x)
                   - ((df1 + df2) / 2.0) * std::log(df1 * x + df2)
                   - lbeta(df1 / 2.0, df2 / 2.0);
 
    return std::exp(log_pdf);
}
 
inline double f_cdf(double x, double df1, double df2)
{
    if (df1 <= 0.0) {
        throw std::invalid_argument("statcpp::f_cdf: df1 must be positive");
    }
    if (df2 <= 0.0) {
        throw std::invalid_argument("statcpp::f_cdf: df2 must be positive");
    }
    if (x <= 0.0) return 0.0;
 
    double z = df1 * x / (df1 * x + df2);
    return betainc(df1 / 2.0, df2 / 2.0, z);
}
 
inline double f_quantile(double p, double df1, double df2)
{
    if (df1 <= 0.0) {
        throw std::invalid_argument("statcpp::f_quantile: df1 must be positive");
    }
    if (df2 <= 0.0) {
        throw std::invalid_argument("statcpp::f_quantile: df2 must be positive");
    }
    if (p < 0.0 || p > 1.0) {
        throw std::invalid_argument("statcpp::f_quantile: p must be in [0, 1]");
    }
    if (p == 0.0) return 0.0;
    if (p == 1.0) return std::numeric_limits<double>::infinity();
 
    // Initial guess using beta quantile
    double z = betaincinv(df1 / 2.0, df2 / 2.0, p);
    double x = df2 * z / (df1 * (1.0 - z));
    if (!std::isfinite(x) || x <= 0.0) {
        x = 1.0;
    }
 
    const double eps = 1e-10;
    const int max_iter = 50;
 
    for (int i = 0; i < max_iter; ++i) {
        if (!std::isfinite(x) || x <= 0.0) {
            x = 1.0;
        }
        double f = f_cdf(x, df1, df2) - p;
        if (std::abs(f) < eps) {
            return x;
        }
        double fprime = f_pdf(x, df1, df2);
        if (fprime == 0.0) break;
 
        double x_new = x - f / fprime;
        if (x_new <= 0.0) {
            x_new = x / 2.0;
        }
 
        if (std::abs(x_new - x) < eps * x) {
            return x_new;
        }
 
        x = x_new;
    }
 
    if (!std::isfinite(x) || x <= 0.0) {
        x = 1.0;
    }
    return x;
}
 
template <typename Engine = default_random_engine>
double f_rand(double df1, double df2, Engine& engine)
{
    if (df1 <= 0.0) {
        throw std::invalid_argument("statcpp::f_rand: df1 must be positive");
    }
    if (df2 <= 0.0) {
        throw std::invalid_argument("statcpp::f_rand: df2 must be positive");
    }
    std::fisher_f_distribution<double> dist(df1, df2);
    return dist(engine);
}
 
inline double f_rand(double df1, double df2)
{
    return f_rand(df1, df2, get_random_engine());
}
 
// ============================================================================
// Log-normal Distribution
// ============================================================================
 
inline double lognormal_pdf(double x, double mu = 0.0, double sigma = 1.0)
{
    if (sigma <= 0.0) {
        throw std::invalid_argument("statcpp::lognormal_pdf: sigma must be positive");
    }
    if (x <= 0.0) {
        return 0.0;
    }
    double log_x = std::log(x);
    double z = (log_x - mu) / sigma;
    return std::exp(-0.5 * z * z) / (x * sigma * sqrt_2_pi);
}
 
inline double lognormal_cdf(double x, double mu = 0.0, double sigma = 1.0)
{
    if (sigma <= 0.0) {
        throw std::invalid_argument("statcpp::lognormal_cdf: sigma must be positive");
    }
    if (x <= 0.0) {
        return 0.0;
    }
    double log_x = std::log(x);
    return norm_cdf((log_x - mu) / sigma);
}
 
inline double lognormal_quantile(double p, double mu = 0.0, double sigma = 1.0)
{
    if (sigma <= 0.0) {
        throw std::invalid_argument("statcpp::lognormal_quantile: sigma must be positive");
    }
    if (p <= 0.0 || p >= 1.0) {
        if (p == 0.0) return 0.0;
        if (p == 1.0) return std::numeric_limits<double>::infinity();
        throw std::invalid_argument("statcpp::lognormal_quantile: p must be in (0, 1)");
    }
    return std::exp(mu + sigma * norm_quantile(p));
}
 
template <typename Engine = default_random_engine>
double lognormal_rand(double mu, double sigma, Engine& engine)
{
    if (sigma <= 0.0) {
        throw std::invalid_argument("statcpp::lognormal_rand: sigma must be positive");
    }
    std::lognormal_distribution<double> dist(mu, sigma);
    return dist(engine);
}
 
inline double lognormal_rand(double mu = 0.0, double sigma = 1.0)
{
    return lognormal_rand(mu, sigma, get_random_engine());
}
 
// ============================================================================
// Weibull Distribution
// ============================================================================
 
inline double weibull_pdf(double x, double shape, double scale = 1.0)
{
    if (shape <= 0.0) {
        throw std::invalid_argument("statcpp::weibull_pdf: shape must be positive");
    }
    if (scale <= 0.0) {
        throw std::invalid_argument("statcpp::weibull_pdf: scale must be positive");
    }
    if (x < 0.0) {
        return 0.0;
    }
    if (x == 0.0) {
        if (shape < 1.0) {
            return std::numeric_limits<double>::infinity();
        } else if (shape == 1.0) {
            return 1.0 / scale;
        } else {
            return 0.0;
        }
    }
    double z = x / scale;
    return (shape / scale) * std::pow(z, shape - 1.0) * std::exp(-std::pow(z, shape));
}
 
inline double weibull_cdf(double x, double shape, double scale = 1.0)
{
    if (shape <= 0.0) {
        throw std::invalid_argument("statcpp::weibull_cdf: shape must be positive");
    }
    if (scale <= 0.0) {
        throw std::invalid_argument("statcpp::weibull_cdf: scale must be positive");
    }
    if (x <= 0.0) {
        return 0.0;
    }
    double z = x / scale;
    return 1.0 - std::exp(-std::pow(z, shape));
}
 
inline double weibull_quantile(double p, double shape, double scale = 1.0)
{
    if (shape <= 0.0) {
        throw std::invalid_argument("statcpp::weibull_quantile: shape must be positive");
    }
    if (scale <= 0.0) {
        throw std::invalid_argument("statcpp::weibull_quantile: scale must be positive");
    }
    if (p < 0.0 || p > 1.0) {
        throw std::invalid_argument("statcpp::weibull_quantile: p must be in [0, 1]");
    }
    if (p == 0.0) return 0.0;
    if (p == 1.0) return std::numeric_limits<double>::infinity();
    return scale * std::pow(-std::log(1.0 - p), 1.0 / shape);
}
 
template <typename Engine = default_random_engine>
double weibull_rand(double shape, double scale, Engine& engine)
{
    if (shape <= 0.0) {
        throw std::invalid_argument("statcpp::weibull_rand: shape must be positive");
    }
    if (scale <= 0.0) {
        throw std::invalid_argument("statcpp::weibull_rand: scale must be positive");
    }
    std::weibull_distribution<double> dist(shape, scale);
    return dist(engine);
}
 
inline double weibull_rand(double shape, double scale = 1.0)
{
    return weibull_rand(shape, scale, get_random_engine());
}
 
// ============================================================================
// Studentized Range Distribution
// ============================================================================
 
inline double studentized_range_cdf(double q, double k, double df)
{
    if (k < 2.0) {
        throw std::invalid_argument("statcpp::studentized_range_cdf: k must be >= 2");
    }
    if (df <= 0.0) {
        throw std::invalid_argument("statcpp::studentized_range_cdf: df must be positive");
    }
    if (q <= 0.0) return 0.0;
 
    // k=2 special case: exact via t-distribution
    // F(q; 2, v) = 2 * t_cdf(q / sqrt(2), v) - 1
    if (k == 2.0) {
        double result = 2.0 * t_cdf(q / sqrt_2, df) - 1.0;
        return std::max(0.0, std::min(1.0, result));
    }
 
    // 12-point Gauss-Legendre half-nodes and half-weights (for inner integral)
    // Following R ptukey.c: only half-nodes stored, symmetric evaluation
    static constexpr int ihalf = 6;
    static constexpr double xleg[ihalf] = {
        0.981560634246719250690549090149,
        0.904117256370474856678465866119,
        0.769902674194304687036893833213,
        0.587317954286617447296702418941,
        0.367831498998180193752691536644,
        0.125233408511468915472441369464
    };
    static constexpr double aleg[ihalf] = {
        0.047175336386511827194615961485,
        0.106939325995318430960254718194,
        0.160078328543346226334652529543,
        0.203167426723065921749064455810,
        0.233492536538354808760849898925,
        0.249147045813402785000562436043
    };
 
    // Inner function: range probability (Hartley's form)
    // Faithful port of R's wprob()
    auto wprob = [&](double w, double cc) -> double {
        static constexpr double C1 = -30.0;
        static constexpr double C2 = -50.0;
        static constexpr double C3 = 60.0;
        static constexpr double bb = 8.0;
        static constexpr double wlar = 3.0;
        static constexpr double wincr1 = 2.0;
        static constexpr double wincr2 = 3.0;
        static constexpr double M_1_SQRT_2PI = 0.398942280401432677939946059934;
 
        double qsqz = w * 0.5;
        if (qsqz >= bb) return 1.0;
 
        // First term: [2*Phi(w/2) - 1]^cc
        double pr_w = 2.0 * norm_cdf(qsqz) - 1.0;
        if (pr_w >= std::exp(C2 / cc))
            pr_w = std::pow(pr_w, cc);
        else
            pr_w = 0.0;
 
        // Interval subdivision for second integral
        double wincr = (w > wlar) ? wincr1 : wincr2;
 
        // Integrate from w/2 to 8 using Gauss-Legendre
        double blb = qsqz;
        double binc = (bb - qsqz) / wincr;
        double bub = blb + binc;
        double einsum = 0.0;
        double cc1 = cc - 1.0;
 
        for (int wi = 0; wi < static_cast<int>(wincr); ++wi) {
            double elsum = 0.0;
            double a = 0.5 * (bub + blb);
            double b = 0.5 * (bub - blb);
 
            for (int jj = 1; jj <= 12; ++jj) {
                double xx;
                int j;
                if (ihalf < jj) {
                    j = 12 - jj;
                    xx = xleg[j];
                } else {
                    j = jj - 1;
                    xx = -xleg[j];
                }
                double ac = a + b * xx;
 
                double qexpo = ac * ac;
                if (qexpo > C3) break;
 
                double pplus = 2.0 * norm_cdf(ac);
                double pminus = 2.0 * norm_cdf(ac - w);
                double rinsum = (pplus * 0.5) - (pminus * 0.5);
 
                if (rinsum >= std::exp(C1 / cc1)) {
                    rinsum = aleg[jj <= ihalf ? jj - 1 : 12 - jj]
                             * std::exp(-0.5 * qexpo) * std::pow(rinsum, cc1);
                    elsum += rinsum;
                }
            }
            elsum *= (2.0 * b) * cc * M_1_SQRT_2PI;
            einsum += elsum;
            blb = bub;
            bub += binc;
        }
 
        pr_w += einsum;
        if (pr_w <= std::exp(C1)) return 0.0;
 
        pr_w = std::pow(pr_w, 1.0);  // rr=1 for Tukey HSD
        return std::max(0.0, std::min(1.0, pr_w));
    };
 
    // Asymptotic case: df > 25000
    if (df > 25000.0) {
        return wprob(q, k);
    }
 
    // General case: integrate over chi-squared density
    // 16-point Gauss-Legendre half-nodes and half-weights (outer integral)
    static constexpr int ihalfq = 8;
    static constexpr double xlegq[ihalfq] = {
        0.989400934991649932596154244360,
        0.944575023073232576077988415535,
        0.865631202387831743880467897712,
        0.755404408355003033895101194847,
        0.617876244402643748446671764049,
        0.458016777657227386342419442984,
        0.281603550779258913230460501460,
        0.095012509837637440185319335425
    };
    static constexpr double alegq[ihalfq] = {
        0.027152459411754094851780572456,
        0.062253523938647892862843836994,
        0.095158511682492784809925107602,
        0.124628971255533872052476282192,
        0.149595988816576732081501730547,
        0.169156519395002538189312079030,
        0.182603415044923588866763667969,
        0.189450610455068496285396723208
    };
 
    // Precompute chi-squared density constants (following R ptukey.c)
    double f2 = df * 0.5;
    double f2lf = (f2 * std::log(df)) - (df * std::log(2.0)) - lgamma(f2);
    double f21 = f2 - 1.0;
    double ff4 = df * 0.25;
 
    // Adaptive interval width
    double ulen;
    if (df <= 100.0) ulen = 1.0;
    else if (df <= 800.0) ulen = 0.5;
    else if (df <= 5000.0) ulen = 0.25;
    else ulen = 0.125;
 
    f2lf += std::log(ulen);
 
    double ans = 0.0;
 
    for (int i = 1; i <= 50; ++i) {
        double otsum = 0.0;
        double twa1 = (2.0 * i - 1.0) * ulen;
 
        for (int jj = 1; jj <= 16; ++jj) {
            double t1;
            int j;
            double t_point;
 
            if (ihalfq < jj) {
                j = jj - ihalfq - 1;
                t_point = twa1 + xlegq[j] * ulen;
                t1 = f2lf + f21 * std::log(t_point) - t_point * ff4;
            } else {
                j = jj - 1;
                t_point = twa1 - xlegq[j] * ulen;
                t1 = f2lf + f21 * std::log(t_point) - t_point * ff4;
            }
 
            if (t_point <= 0.0) continue;
            if (t1 >= -30.0) {
                double qsqz;
                if (ihalfq < jj) {
                    qsqz = q * std::sqrt((xlegq[j] * ulen + twa1) * 0.5);
                } else {
                    qsqz = q * std::sqrt((-(xlegq[j] * ulen) + twa1) * 0.5);
                }
 
                double wprb = wprob(qsqz, k);
                otsum += (wprb * alegq[j]) * std::exp(t1);
            }
        }
 
        // Convergence check (R: if i*ulen >= 1.0 && otsum <= 1e-14)
        if (static_cast<double>(i) * ulen >= 1.0 && otsum <= 1e-14) {
            break;
        }
 
        ans += otsum;
    }
 
    return std::max(0.0, std::min(1.0, ans));
}
 
inline double studentized_range_quantile(double p, double k, double df)
{
    if (k < 2.0) {
        throw std::invalid_argument("statcpp::studentized_range_quantile: k must be >= 2");
    }
    if (df <= 0.0) {
        throw std::invalid_argument("statcpp::studentized_range_quantile: df must be positive");
    }
    if (p < 0.0 || p > 1.0) {
        throw std::invalid_argument("statcpp::studentized_range_quantile: p must be in [0, 1]");
    }
    if (p == 0.0) return 0.0;
    if (p == 1.0) return std::numeric_limits<double>::infinity();
 
    // Initial guess: heuristic based on normal quantile
    double q0 = norm_quantile(1.0 - (1.0 - p) * 0.5);
    double x = q0 * (1.0 + 0.2 * (k - 2.0) / std::sqrt(df));
    if (!std::isfinite(x) || x < 1.0) x = 1.0;
 
    // Newton-Raphson iteration
    const double eps = 1e-10;
    const int max_iter = 50;
 
    for (int i = 0; i < max_iter; ++i) {
        if (!std::isfinite(x) || x <= 0.0) {
            x = 1.0;
        }
        double f = studentized_range_cdf(x, k, df) - p;
        if (std::abs(f) < eps) return x;
 
        // Numerical derivative (central difference)
        double dx = std::max(1e-6, x * 1e-6);
        double fprime = (studentized_range_cdf(x + dx, k, df)
                       - studentized_range_cdf(x - dx, k, df)) / (2.0 * dx);
        if (fprime <= 0.0) break;
 
        double x_new = x - f / fprime;
        if (x_new <= 0.0) x_new = x * 0.5;
 
        if (std::abs(x_new - x) < eps * (1.0 + std::abs(x))) return x_new;
        x = x_new;
    }
 
    if (!std::isfinite(x) || x <= 0.0) {
        x = 1.0;
    }
    return x;
}
 
} // namespace statcpp