t_distribution_cdf

C++ Function Reference

1 Signature

double BigDataStatMeth::t_distribution_cdf(double t, double df)

2 Description

Compute p-value for correlation coefficient using t-distribution.

3 Parameters

• r (``): Correlation coefficient [-1, 1]
• n (``): Sample size (number of observations)
• t (double): T-statistic value
• df (double): Degrees of freedom (ν)

4 Returns

Two-tailed p-value [0, 1], or NaN if invalid input

5 Details

Computes the two-tailed p-value for a correlation coefficient using the t-statistic: t = r * sqrt((n-2)/(1-r²))

6 Caller Graph

7 Source Code

NoteImplementation

File: inst/include/hdf5Algebra/matrixCorrelation.hpp • Lines 236-338

inline double t_distribution_cdf(double t, double df) {
        if (df <= 0) return std::numeric_limits<double>::quiet_NaN();
        if (std::isinf(t)) return (t > 0) ? 1.0 : 0.0;
        if (t == 0) return 0.5;
        
        // For large df, use normal approximation
        if (df > 100) {
            // Standard normal CDF approximation
            double z = t;
            return 0.5 * (1.0 + std::erf(z / std::sqrt(2.0)));
        }
        
        // Use the relationship: t-CDF = 0.5 + sign(t) * I_x(1/2, df/2) / 2
        // where x = t²/(t² + df) and I is the regularized incomplete beta function
        double x = t * t / (t * t + df);
        
        // Compute incomplete beta function I_x(1/2, df/2) using continued fraction
        // This is equivalent to R's implementation
        double a = 0.5;
        double b = df / 2.0;
        
        // For small x, use series expansion
        if (x < (a + 1.0) / (a + b + 2.0)) {
            // Use continued fraction for B(x; a, b)
            double bt = std::exp(std::lgamma(a + b) - std::lgamma(a) - std::lgamma(b) + 
                                 a * std::log(x) + b * std::log(1.0 - x));
            
            // Continued fraction approximation
            double qab = a + b;
            double qap = a + 1.0;
            double qam = a - 1.0;
            double c = 1.0;
            double d = 1.0 - qab * x / qap;
            
            if (std::abs(d) < 1e-30) d = 1e-30;
            d = 1.0 / d;
            double h = d;
            
            for (int m = 1; m <= 100; ++m) {
                int m2 = 2 * m;
                double aa = m * (b - m) * x / ((qam + m2) * (a + m2));
                d = 1.0 + aa * d;
                if (std::abs(d) < 1e-30) d = 1e-30;
                c = 1.0 + aa / c;
                if (std::abs(c) < 1e-30) c = 1e-30;
                d = 1.0 / d;
                h *= d * c;
                
                aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
                d = 1.0 + aa * d;
                if (std::abs(d) < 1e-30) d = 1e-30;
                c = 1.0 + aa / c;
                if (std::abs(c) < 1e-30) c = 1e-30;
                d = 1.0 / d;
                double del = d * c;
                h *= del;
                
                if (std::abs(del - 1.0) < 1e-12) break;
            }
            
            double beta_incomplete = bt * h / a;
            return 0.5 + (t > 0 ? 1 : -1) * beta_incomplete / 2.0;
            
        } else {
            // Use the symmetry relation
            double bt = std::exp(std::lgamma(a + b) - std::lgamma(a) - std::lgamma(b) + 
                                 b * std::log(1.0 - x) + a * std::log(x));
            
            // Continued fraction for the complementary case
            double qab = a + b;
            double c = 1.0;
            double d = 1.0 - qab * (1.0 - x) / (b + 1.0);
            
            if (std::abs(d) < 1e-30) d = 1e-30;
            d = 1.0 / d;
            double h = d;
            
            for (int m = 1; m <= 100; ++m) {
                int m2 = 2 * m;
                double aa = m * (a - m) * (1.0 - x) / ((b - 1.0 + m2) * (b + m2));
                d = 1.0 + aa * d;
                if (std::abs(d) < 1e-30) d = 1e-30;
                c = 1.0 + aa / c;
                if (std::abs(c) < 1e-30) c = 1e-30;
                d = 1.0 / d;
                h *= d * c;
                
                aa = -(b + m) * (qab + m) * (1.0 - x) / ((b + m2) * (b + 1.0 + m2));
                d = 1.0 + aa * d;
                if (std::abs(d) < 1e-30) d = 1e-30;
                c = 1.0 + aa / c;
                if (std::abs(c) < 1e-30) c = 1e-30;
                d = 1.0 / d;
                double del = d * c;
                h *= del;
                
                if (std::abs(del - 1.0) < 1e-12) break;
            }
            
            double beta_incomplete = 1.0 - bt * h / b;
            return 0.5 + (t > 0 ? 1 : -1) * beta_incomplete / 2.0;
        }
    }

8 Usage Example

#include "BigDataStatMeth.hpp"

// Example usage
auto result = t_distribution_cdf(...);