1.   定理

对于实数 $a,b>0$ 有贝塔函数定义：

Gamma 函数由 Leonhard Euler 于 1730 年“内插阶乘”引入，定义为：

其中 $\operatorname{Re}(x)>0$ 保证积分收敛，被积函数 $\mathrm{e}^{-t}t^{x-1}$ 在 $t\to0^{+}$ 与 $t\to+\infty$ 两端均趋于零。

Weierstrass 无穷乘积是 Gamma 函数的倒数，被定义为：

其中 $\gamma\approx0.5772$ 为 Euler-Mascheroni 常数，由此 $\dfrac{1}{\Gamma(x)}$ 是整函数，$\Gamma(x)$ 在 $x=0,-1,-2,\ldots$ 处有一阶极点，无零点。

分部积分定义：

洛必达法则有：

若 $x\in\mathbb{N}$ 则 Gamma 数列被定义为：

2.   正则化不完全贝塔函数

将 $\beta$ 函数的积分上限换成变量 $x$，其中 $0\leqslant x\leqslant 1$，得到不完全贝塔函数：

概率统计中，需保证总概率恒为 1，因此将不完全 $\beta$ 函数除以 $\beta$ 函数进行归一化，其中 $a,b\in\mathbb{N}$，得到正则化不完全 $\beta$ 函数：

3.   计算对数伽马

Gamma 函数增长极快，第一种情况是利用 $\Gamma(x+1)=x\Gamma(x)$ 递推，将伽马变到 $1\sim2$ 范围，但比如 $\Gamma(172)=171!>10^{308}$，即超过 double 双精度浮点数最大范围。第二种是考虑到使用对数压缩范围，这样可以使得计算不超过浮点数范围。

对于正数，使用斯特林近似：

对于负数，使用欧拉反射公式：

4.   无穷连分式展开

根据美国国家标准与技术研究院 NIST 出版的数学函数数字图书库 Digital Library of Mathematical Functions 权威定义，正则化不完全 $\beta$ 函数可以通过以下无限连分数展开进行超高精度逼近：

偶数项系数序列：

奇数项系数序列：

连分数收敛条件：

评估点 $x$ 满足此核心不等式时，概率积分的几何权重分布主要集中在评估区间的前半部分，连分数的逼近截断误差呈现出极其优异的指数级衰减。然而，当输入的 $x$ 不满足条件时，连分数分母中的相互抵消效应会严重恶化，导致收敛速度急剧下降，如果强行继续迭代，不仅会消耗大量的 CPU 计算时钟周期，甚至可能由于双精度浮点数 double 类型尾数位 Mantissa 的截断误差累积，导致最终计算彻底发散失效。

为从根本上解决这一计算理论瓶颈，利用正则化不完全 $\beta$ 函数的镜像对称性等价恒等式：

任意形状的 $\beta$ 分布曲线，位于左侧尾部的概率累积面积，严格等价于其形状参数相互对调后，位于右侧相应镜像区间的尾部概率面积之补集。

实际编程算法实现层面，当检测到输入变量 $x$ 超出收敛域时，程序会立刻自动触发参数翻转机制，将目标计算任务转换为计算 $I(1-x;b,a)$，由于 $x$ 与 $1-x$ 处关于区间中点的互补状态，这种基于镜像对称性的动态条件分支策略，达成一种优雅的自适应保护机制，确保无论外部调用栈传入的输入参数如何畸形或极端，底层最深处的连分数引擎始终且必须运行在性能最佳的指数级收敛区间内。

5.   函数序列的迭代 Lentz 算法

Lentz 算法是一种高效且数值稳定计算连分数的方法，由 William J. Lentz 在 1976 年，用于计算贝塞尔函数等特殊函数的值，但后来被广泛用于各种无穷连分数的数值计算中，Lentz 算法求解的目标是形如：

连分数直接从上到下计算容易出现除零或者数值溢出问题，Lentz 算法通过迭代更新避免了这些问题。Lentz 算法核心思想：

• 使用辅助状态 $C_n$ 和 $D_n$ 迭代更新分子分母的比值；
• 引入一个很小的数 $\varepsilon$ 防止除零；
• 每次迭代更新当前的近似值 $F_n$ 直到收敛。

计算正则化不完全 $\beta$ 函数，其中 $b_0=b_1=b_2=\cdots=1$，而 $a_i$ 严格对应 $d_i$，并且设置 Lentz 初始值：

6.   高斯积分

记 $A=\displaystyle\int_{-\infty}^{\infty}\mathrm{e}^{\displaystyle-x^2}\mathrm{d}x$，则：

令：

根据雅可比行列式可得：

则：

记 $u=r^2$ 则微分：

则：

根据高斯积分，易算得有如下推论：

7.   计算 Student $t$ 分布的累积分布函数

Student $t$ 的 PDF 概率密度函数：

由于：

则有：

则有累积概率分布函数：

令 $u=\dfrac{\nu}{\nu+x^2}$，则有：

当 $x\to\infty$ 时，$u\to0$，当 $x\to t$ 时，$u\to\dfrac{\nu}{\nu+t^2}$，记 $s=\dfrac{\nu}{\nu+t^2}$，则有：

因此 Student $t$ 分布的 CDF 结果如下：

自由度 ν = 1

自由度 ν = 5

自由度 ν = 30

8.   全局平滑映射

为了消除编程中 if-else 分支判断正负，加快 CPU 计算，我们把 $t\in(-\infty,+\infty)$ 映射到正则化不完全 $\beta$ 函数所需的 $x\in[0,1]$ 区间，刚才是令 $s=\dfrac{\nu}{\nu+t^2}$，但代码里必须进行 if t > 0 和 if t <= 0 分段处理，并且参数是固定的 $\dfrac{\nu}{2}$ 和 $\dfrac{1}{2}$，现在用另一种映射变量：

这是一种单调递增的映射，有这些特点：

• 当 $t\to-\infty$ 时 $z\to0$；
• 当 $t=0$ 时 $z=0.5$；
• 当 $t\to+\infty$ 时 $z\to1$。

它把整个实数轴平滑地拉伸并压缩到了 $[0,1]$ 区间内，不需要任何 if-else 分支来判断正负。

当 Student $t$ 分布以 $t=0$ 为中心对称时，让正则化不完全 $\beta$ 函数的两个参数都为 $\dfrac{\nu}{2}$，则利用贝塔函数对称性恒等式：

若将 $z=\dfrac{1}{2}\left(1+\dfrac{t}{\sqrt{t^2+\nu}}\right)$ 代入 $1-(1-2z)^2$ 可得：

这与上述 Student $t$ 分布 CDF 在数学推导上是等价的。

9.   C 语言构建 t 分布临界置信水平表

现在有了第 1～8 部分的数学理论，我们容易迭代求得 $t$ 分布的 CDF 累积概率分布，再根据定义可知，$t$ 分布的 PDF 概率密度函数是 CDF 的导数，所以：

为了提高逼近的精度，工程上通常不采用单边差分，而是采用中心差分：

你对比一下，离谱吧，只需要调用两次 incbeta 正则化不完全贝塔函数，再将 $h$ 步长设置极小如 1e-7 浮点数，完全不需要用 $\Gamma$ 函数和复杂的代数解析式，就能得到精确的 PDF 值！

此时，为了构建 $t$ 分布临界值表，我们只需要澄清一个道理：由显著性水平 $\alpha$ 可以得到目标累积概率 $p$，除此之外，咱们还有自由度 $\nu$，如何求未知数 $x$？即求解非线性方程：

考虑 Newton-Raphson 迭代算法：

最后编译代码执行程序得到 Student $t$ 分布临界置信水平表：

C 语言代码运行结果

源码如下：

#include <math.h>
#include <stdio.h>

double incbeta(double x, double a, double b) {
    static const double epsilon = 1.0e-30;
    static const double condition = 1.0e-12;

    if (x < 0.0 || x > 1.0) return INFINITY;
    if (x > (a + 1.0) / (a + b + 2.0)) return 1.0 - incbeta(1.0 - x, b, a);

    double c = 1.0;
    double d = 0.0;
    double f = 1.0;

    const double lbeta_ab = lgamma(a) + lgamma(b) - lgamma(a + b);
    const double front = exp(a * log(x) + b * log(1.0 - x) - lbeta_ab) / a;

    for (int i = 0; i <= 200; ++i) {
        int m = i / 2;
        double numerator;

        if (i == 0) numerator = 1.0;
        else if (i % 2 == 0) numerator = (m * (b - m) * x) / ((a + 2.0 * m - 1.0) * (a + 2.0 * m));
        else numerator = - ((a + m) * (a + b + m) * x) / ((a + 2.0 * m) * (a + 2.0 * m + 1));

        d = 1.0 + numerator * d;
        if (fabs(d) < epsilon) d = epsilon;
        d = 1.0 / d;

        c = 1.0 + numerator / c;
        if (fabs(c) < epsilon) c = epsilon;

        double cd = c * d;
        f = f * cd;
        if (fabs(cd - 1.0) < condition) return front * (f - 1.0);
    }
    return INFINITY;
}

double cdf_student_t(double t, double df) {
    if (df <= 0.0) return NAN;
    if (t == 0.0) return 0.5;
    double cache = sqrt(t * t + df);
    return incbeta((t + cache) / (2.0 * cache), df / 2.0, df / 2.0);
}

double pdf_student_t(double t, double df) {
    double h = 1e-7;
    double cdf_plus = cdf_student_t(t + h, df);
    double cdf_minus = cdf_student_t(t - h, df);
    return (cdf_plus - cdf_minus) / (2.0 * h);
}

double newton_raphson_t(double p, double df) {
    if (p <= 0.0) return -INFINITY;
    if (p >= 1.0) return INFINITY;
    if (p == 0.5) return 0.0;

    double x = 0.0;
    int epochs = 50;
    double epsilon = 1e-12;

    for (int n = 0; n < epochs; n++) {
        double g = cdf_student_t(x, df) - p;
        if (fabs(g) < epsilon) break;

        double g_derivative = pdf_student_t(x, df);
        if (g_derivative == 0.0) break;

        x = x - g / g_derivative;
    }
    return x;
}

void build_student_t_critical_table() {
    double single_alpha[] = {0.25, 0.2, 0.15, 0.1, 0.05, 0.025, 0.01, 0.005, 0.001};
    double size = sizeof(single_alpha) / sizeof(single_alpha[0]);

    printf("===============================================================================================================\n");
    printf("                                        P{t(v) > t(v;1-alpha)} = alpha\n");
    printf("---------------------------------------------------------------------------------------------------------------\n");
    printf(" alpha / 2 |");
    for (int i = 0; i < size; i++) printf(" %7g  |", single_alpha[i]);
    printf("\n");
    printf("   alpha   |");
    for (int i = 0; i < size; i++) printf(" %7g  |", single_alpha[i] * 2.0);
    printf("\n");
    printf("---------------------------------------------------------------------------------------------------------------\n");
    printf("   df(v)   |\n");
    printf("---------------------------------------------------------------------------------------------------------------\n");
    for (double v = 1; v < 16; v++) {
        printf("    %2.0lf     |", v);
        for (int i = 0; i < size; i++) {
            double p = 1.0 - single_alpha[i];
            double t = newton_raphson_t(p, v);
            printf(" %8.4lf |", t);
        }
        printf("\n");
    }
    printf("===============================================================================================================\n");
}

int main(int argc, char *argv[], char *envs[]) {
    build_student_t_critical_table();
    return 0;
}

彩蛋：补充

F 分布的累积分布函数：

其中：

同理类推 Student $t$ 分布计算临界值方法，可以快速求出 F 分布的临界值表。

double cdf_f(double F, double df1, double df2) {
    if (F <= 0) return 0.0;
    double x = (df1 * F) / (df1 * F + df2);
    return incbeta(x, df1 / 2.0, df2 / 2.0);
}

double pdf_f(double F, double df1, double df2) {
    double h = 1e-7;
    double cdf_minus = (F - h < 0.0) ? 0.0 : cdf_f(F - h, df1, df2);
    return (cdf_f(F + h, df1, df2) - cdf_minus) / (2.0 * h);
}

double newton_raphson_f(double p, double df1, double df2) {
    if (p <= 0.0) return 0.0;
    if (p >= 1.0) return INFINITY;

    double x = 1.0;
    int epochs = 50;
    double epsilon = 1e-12;

    for (int n = 0; n < epochs; n++) {
        double g = cdf_f(x, df1, df2) - p;
        if (fabs(g) < epsilon) break;

        double g_derivative = pdf_f(x, df1, df2);
        if (g_derivative == 0.0) break;

        x = x - g / g_derivative;
        if (x <= 0.0) x = 0.0000001;
    }
    return x;
}

卡方 $\chi^2$ 分布衡量是多个标准正态分布的平方和，概率密度函数定义为：

当 F 分布自由度 $\nu_2\to\infty$ 时，可得 $\chi^2$ 关系恒等式：

double cdf_approx_chi_square(double p, double df) {
    double infinity_df2 = 1e6;
    return df * newton_raphson_f(p, df, infinity_df2);
}

测试结果：

int main(int argc, char *argv[], char *envs[]) {
    double alpha = 0.01;
    double p = 1.0 - alpha;
    double df1 = 5;
    double df2 = 10;
    double df = 7;

    printf("alpha = %.3lf, df1 = %.0f, df2 = %.0f, F = %.4lf\n", alpha, df1, df2, newton_raphson_f(p, df1, df2));
    printf("alpha = %.3lf, df = %.0f, chi^2 = %.4lf\n", alpha, df, cdf_approx_chi_square(p, df));
    return 0;
}

alpha = 0.010, df1 = 5, df2 = 10, F = 5.6363

alpha = 0.010, df = 7, chi^2 = 18.4754