C++：评估非中心 T 的累积密度函数（CDF）概率密度函数（PDF）（附带源码）

一、项目背景详细介绍

在统计推断、信号处理、医学统计、金融计量、A/B 测试以及假设检验中，Student’s t 分布是一个极其重要的概率分布。
当总体方差未知、样本量有限时，t 分布是正态分布的重要替代。

然而，在很多真实工程与科研场景中，零假设并不总是“均值为 0”。
例如：

信号检测中存在已知偏移量
医学实验中药物效应不为零
金融中收益率存在漂移
功效分析（Power Analysis）
贝叶斯推断中的后验预测

这时，**非中心 t 分布（Non-central t Distribution）**就成为核心工具。

1.1 什么是非中心 t 分布

1.2 为什么 C++ 中“很难”直接算非中心 t？

现实问题：

标准库<cmath>没有非中心 t
Boost 虽然有实现，但：
- 有时不能用第三方库
- 有时需要“透明可控”的数值实现
- 有时用于教学/科研验证

因此，手写一个数值稳定、可读、可扩展的实现具有很高价值。

二、项目需求详细介绍

本项目目标如下：

2.1 功能需求

实现一个C++ 数值库级别模块，支持：

非中心 t 分布 PDF
非中心 t 分布 CDF
支持任意：

2.2 数值需求

稳定性高（避免下溢/上溢）
对中等自由度（1 ~ 1000）表现良好
可控截断误差
不依赖第三方库

三、相关技术详细介绍

3.1 非中心 t 的数学定义

PDF（概率密度函数）

3.2 中心 t 分布基础

3.3 数值计算核心技巧

技术	用途
`lgamma`	避免 Gamma 溢出
Poisson 截断	有限级数逼近
对称性	提升稳定性
log-domain	防止下溢

四、实现思路详细介绍

4.1 总体架构

noncentral_t/ ├── math_utils │ ├── gamma │ ├── beta │ └── poisson ├── student_t │ ├── central_pdf │ ├── central_cdf │ ├── noncentral_pdf │ └── noncentral_cdf └── main

4.2 非中心 t PDF 实现思路

计算 Poisson 权重
对每个 k：
- 计算自由度 ν+2k\nu + 2kν+2k
- 调用中心 t PDF
累加直到：
- 权重足够小
- 或达到最大迭代次数

4.3 非中心 t CDF 实现思路

与 PDF 类似：

核心是中心 t CDF
利用级数展开
对称处理负值

五、完整实现代码

/***************************************************** * File: noncentral_t.cpp * Description: * Evaluation of non-central Student's t * PDF and CDF in pure C++ * Standard: C++17 *****************************************************/ #include <cmath> #include <iostream> #include <limits> namespace math { /*************** Gamma utilities ****************/ double log_gamma(double x) { return std::lgamma(x); } /*************** Beta function ******************/ double incomplete_beta(double a, double b, double x) { // 简化实现：数值积分（教学用途） const int N = 2000; double h = x / N; double sum = 0.0; for (int i = 1; i < N; ++i) { double t = i * h; sum += std::pow(t, a - 1) * std::pow(1 - t, b - 1); } return sum * h; } double regularized_beta(double a, double b, double x) { double beta = std::exp(log_gamma(a) + log_gamma(b) - log_gamma(a + b)); return incomplete_beta(a, b, x) / beta; } /*************** Central t PDF ******************/ double central_t_pdf(double t, double nu) { double logc = log_gamma((nu + 1) / 2.0) - log_gamma(nu / 2.0) - 0.5 * std::log(nu * M_PI); double logp = logc - (nu + 1) / 2.0 * std::log(1 + t * t / nu); return std::exp(logp); } /*************** Central t CDF ******************/ double central_t_cdf(double t, double nu) { if (t == 0.0) return 0.5; double x = nu / (nu + t * t); double ib = regularized_beta(nu / 2.0, 0.5, x); if (t > 0) return 1.0 - 0.5 * ib; else return 0.5 * ib; } /*************** Poisson weight *****************/ double poisson_weight(int k, double lambda) { return std::exp(k * std::log(lambda) - lambda - log_gamma(k + 1)); } } // namespace math /*************** Non-central t PDF ***************/ double noncentral_t_pdf(double t, double nu, double delta) { const int MAX_K = 100; double lambda = delta * delta / 2.0; double sum = 0.0; for (int k = 0; k < MAX_K; ++k) { double w = math::poisson_weight(k, lambda); double df = nu + 2 * k; sum += w * math::central_t_pdf(t, df); if (w < 1e-12) break; } return sum; } /*************** Non-central t CDF ***************/ double noncentral_t_cdf(double t, double nu, double delta) { const int MAX_K = 100; double lambda = delta * delta / 2.0; double sum = 0.0; for (int k = 0; k < MAX_K; ++k) { double w = math::poisson_weight(k, lambda); double df = nu + 2 * k; sum += w * math::central_t_cdf(t, df); if (w < 1e-12) break; } return sum; } /*********************** Main ********************/ int main() { double t = 1.5; double nu = 10.0; double delta = 2.0; std::cout << "Non-central t PDF: " << noncentral_t_pdf(t, nu, delta) << std::endl; std::cout << "Non-central t CDF: " << noncentral_t_cdf(t, nu, delta) << std::endl; return 0; }