# Math in Biometrics

# Mathematical Analysis

# Leibniz's Notation

<h3>Who is Leibniz</h3>
<p>Gottfried Wilhelm Leibniz (1646 - 1716) was a 17th century German mathematician. He’s often credited with developing many of the main principles of differential and integral calculus, and is primarily recognized for what we now call Leibniz’s notation.</p>
<div align="center">
  <img src="https://macrolib.com/uploads/images/gallery/2025-11/leibniz600.png">
</div>
<br>

<h3>Leibniz's Notation System</h3>
<p>The derivative of a function based on today’s standard  is given by:</p>
<div align="center">
  <img src="https://latex.codecogs.com/svg.image?\displaystyle\lim_{\Delta&space;x\to&space;0}\frac{\Delta&space;y}{\Delta&space;x}=\displaystyle\lim_{\Delta&space;x\to&space;0}\frac{f(x&plus;\Delta&space;x)-f(x)}{\Delta&space;x}">
</div>
<br>
<p>Leibniz's notation expresses the derivative as:</p>
<div align="center">
  <img src="https://latex.codecogs.com/svg.image?\frac{dy}{dx}=\frac{f(x&plus;dx)-f(x)}{dx}">
</div>
<p>where goes <img src="https://latex.codecogs.com/svg.image?dx" style="vertical-align:middle"> toward 0.</p>

<h3>Fractional Behavior</h3>
<p>Let’s review some examples where Leibniz’s notation is often utilized. The Chain Rule using Lagrange Notation states:</p>
<div align="center">
  <img src="https://latex.codecogs.com/svg.image?h^{'}(x)=f^{'}(g(x))g^{'}(x)">
</div>
<br>
<p>We can translate the above Chain Rule into Leibniz's Notation as:</p>
<div align="center">
  <img src="https://latex.codecogs.com/svg.image?\frac{dy}{dx}=\frac{dy}{du}\frac{du}{dx}">
</div>
<br>
<p>In the above equations, we can see how Leibniz’s Notation behaves similarly to a fraction, although it must be emphasized that the derivative is not a fraction.</p>

<h3>Try It on a Function</h3>
<p>Let us try <img src="https://latex.codecogs.com/svg.image?\inline&space;f(x)=x^{2}" style="vertical-align:middle">.</p>
<details>
<summary>Click to Expand/Collapse</summary>
  <p>Expand the function:</p>
  <div align="center">
    <img src="https://latex.codecogs.com/svg.image?\frac{dy}{dx}=\frac{(x&plus;dx)^{2}-x^{2}}{dx}=\frac{x^{2}&plus;2x(dx)&plus;(dx)^{2}-x^{2}}{dx}">
  </div>
  <p>Simplify fraction:</p>
  <div align="center">
    <img src="https://latex.codecogs.com/svg.image?\frac{dy}{dx}=\frac{2x(dx)&plus;(dx)^{2}}{dx}=2x&plus;dx=2x">
  </div>
<br>
</details>
<p>So the derivative of <img src="https://latex.codecogs.com/svg.image?\inline&space;&space;x^{2}"> is <img src="https://latex.codecogs.com/svg.image?\inline&space;&space;2x">.</p>

# Numerical Computation

# Owen's T/Q 函数在生物统计学样本量计算中的应用

## 1. 使用场景：它解决什么问题

**一句话直觉**：当你要回答"做多少次实验才够，才能有把握地宣称两个东西**几乎一样**（而不是**不一样**）"，普通的"差异检验"功效公式就会失效，这时就需要 Owen's T/Q 函数来精确计算把握度（power）和样本量。

最典型、也最规范的应用场景是**药物生物等效性（bioequivalence, BE）试验**：

- 仿制药（generic drug）上市前，必须证明它与原研药（reference）在吸收速度和程度上**等效**（通常是药代动力学参数 AUC、C_max 的几何均值比落在 $0.80$–$1.25$ 之间）。
- 法规上要求采用 **双单侧检验（TOST）** 程序：把"等效"拆成两个单侧检验，**两个检验都必须显著**才算等效。

**为什么普通功效概念不够用**：在"差异检验"里，功效是"真实差异存在时能检出的概率"，随样本量增大而单调趋近 1；而"等效检验"的功效是"真实比值落在等效区间内时、置信区间能被完整装进 $[0.80,1.25]$ 的概率"，它由**两个共享同一随机分母的非中心 t 统计量**共同决定。这个**联合概率**无法用单一非中心 t 分布直接写出，而恰好可以用 **Owen 的 Q 函数**（以及更基础的 Owen 的 T 函数）来精确表达。

> **Owen's T 与 Owen's Q 的分工**：
> - **Owen's T**（Owen 1956）解决的是**二元正态**矩形区域的概率；
> - **Owen's Q**（Owen 1965）把同样的思想推广到**二元非中心 t 分布**——而 TOST 的功效正是"相关系数为 $-1$ 的二元非中心 t"这种特殊情形的概率（Owen 那篇论文的标题就叫 *A Special Case of a Bivariate Non-Central t-Distribution*）。

**一个关键结论（后文会推导并验证）**：TOST 的精确功效

$$
1-\beta \;=\; F_\nu\!\left(-t_{1-\alpha,\nu};\delta_2\right) - F_\nu\!\left(t_{1-\alpha,\nu};\delta_1\right) - \big[Q\text{-校正项}\big],
$$

其中 $F_\nu(\cdot;\delta)$ 是非中心 t 的 CDF（它本身就是 $Q_\nu$ 在积分下限为 $0$ 时的特例），方括号里是 Owen's Q 在积分下限 $R>0$ 处的小校正。**实用上**这个校正项通常极小，但"必须同时通过两个单侧检验"这个两尾结构带来的效应却很大——这正是必须用 Owen 的 Q 函数而非简单正态近似的根本原因。

---

## 2. 基础知识

### 2.1 非中心 t 分布

**直觉**：普通的（中心）t 分布描述"样本均值对标准误的比值"在**零假设为真**时的行为；非中心 t 分布描述**备择假设为真**（真实效应 $\delta\neq 0$）时同一个比值的行为。功效分析的本质，就是在备择假设下计算这个比值的分布。

**定义**：若 $Z\sim N(0,1)$ 与 $V\sim\chi^2_\nu$ 独立，则

$$
T \;=\; \frac{Z+\delta}{\sqrt{V/\nu}} \;\sim\; t_\nu(\delta),
$$

$\delta\in\mathbb R$ 为**非中心参数（noncentrality parameter, NCP）**，$\nu$ 为自由度。$\delta=0$ 时退化为中心 t 分布。

**CDF 的积分表示**（这是通往 Owen's Q 的关键一步）。设 $W=\sqrt{V/\nu}$，则

$$
F_\nu(t;\delta)=\Pr(T\le t)=\int_0^\infty \Phi(t w-\delta)\,f_W(w)\,dw,
\qquad
f_W(w)=\frac{2(\nu/2)^{\nu/2}}{\Gamma(\nu/2)}\,w^{\nu-1}e^{-\nu w^2/2},\quad w>0.
$$

直观解释：先固定分母 $W=w$，此时 $T$ 条件地服从 $N(\delta/w,\,1/w^2)$ 形态，$\Pr(T\le t\mid w)=\Phi(tw-\delta)$；再对 $w$ 的（卡方派生）密度积分。

**矩**（对 $\nu>1$、$\nu>2$ 分别成立）：

$$
\mathbb E[T]=\delta\sqrt{\frac{\nu}{2}}\,\frac{\Gamma\!\big(\frac{\nu-1}{2}\big)}{\Gamma\!\big(\frac{\nu}{2}\big)},
\qquad
\operatorname{Var}(T)=\frac{\nu}{\nu-2}(1+\delta^2)-\big(\mathbb E[T]\big)^2.
$$

注意 $\mathbb E[T]\approx\delta$（当 $\nu$ 较大时），这解释了为什么 NCP $\delta$ 可以近似看成"检验统计量期望偏离零假设的量"。

### 2.2 Owen's T 函数（二元正态概率）

**直觉**：要算"两个都服从正态的变量 $X,Y$ **同时**不超过各自阈值"的概率（二元正态矩形概率），不能简单相乘（它们相关），Owen's T 就是那个"修正相关性的积分"。

**定义**（Owen 1956）：

$$
T(h,a)=\frac{1}{2\pi}\int_0^{a}\frac{\exp\!\big[-\tfrac{h^2}{2}(1+x^2)\big]}{1+x^2}\,dx,
\qquad -\infty<h<\infty,\ -\infty<a<\infty.
$$

**基本性质**：

$$
T(h,0)=0,\qquad
T(0,a)=\frac{1}{2\pi}\arctan a,\qquad
T(-h,a)=T(h,a),\qquad
T(h,-a)=-T(h,a),\qquad
T(h,\infty)=\tfrac12\,\Phi(-|h|).
$$

**与二元正态的联系**：对相关系数为 $\rho$ 的标准二元正态，

$$
\Pr(X\le h,\ Y\le k)=\Phi(h)\Phi(k)+\frac{1}{2\pi}\int_0^{\rho}\frac{1}{\sqrt{1-r^2}}\exp\!\left[-\frac{h^2-2rhk+k^2}{2(1-r^2)}\right]dr,
$$

Owen (1956) 进一步把它约化为 T 函数：

$$
\Pr(X\le h,Y\le k)=\tfrac12\Phi(h)+\tfrac12\Phi(k)-T\!\left(h,\frac{k-\rho h}{h\sqrt{1-\rho^2}}\right)-T\!\left(k,\frac{h-\rho k}{k\sqrt{1-\rho^2}}\right)-c,
\qquad c=\begin{cases}0,&hk\ge0\\ \tfrac12,&hk<0.\end{cases}
$$

（$h$ 或 $k=0$ 为退化边界，需单独处理。）

### 2.3 Owen's Q 函数（二元非中心 t / 非中心 t CDF）

**直觉**：把 2.2 的"二元正态矩形概率"换成"二元非中心 t"，并且考虑两个统计量**共享同一个分母**（相关系数 $\pm 1$）这种特殊情形，得到的推广就是 Owen's Q。

**定义**（Owen 1965；本文采用"约定 A"，即积分核中为 $tx/\nu$）：

$$
Q_\nu(t,\delta;0,R)=\frac{\sqrt{2\pi}}{\Gamma\!\big(\tfrac\nu2\big)\,2^{(\nu-2)/2}}\int_R^{\infty}\Phi\!\left(\frac{t x}{\nu}-\delta\right)x^{\nu-1}\phi(x)\,dx,
\qquad \phi(x)=\frac{e^{-x^2/2}}{\sqrt{2\pi}}.
$$

（不同文献可能把核写成 $tx/\sqrt\nu$，二者只差一个换元 $x\mapsto x\sqrt\nu$，读文献时务必核对约定。）

**关键恒等式**：作换元 $x=w\sqrt\nu$ 后，$Q_\nu$ 与 2.1 的非中心 t CDF 一致，即

$$
F_\nu(t;\delta)=Q_\nu\!\big(t\sqrt\nu,\ \delta;\ 0,\ 0\big).
$$

也就是说：**非中心 t 的 CDF 就是 Owen's Q 函数在积分下限 $R=0$ 时的特例**；而 $R>0$ 时的 $Q_\nu$ 用来处理"截断积分"，正是 TOST 功效里需要的那个校正。

### 2.4 生物等效性的假设检验（TOST）

**背景与尺度**：药代参数通常近似对数正态，故分析在**对数尺度**进行。记

$$
\mu=\ln\!\left(\frac{\text{test}}{\text{reference}}\right)=\ln(\text{真实比值}),
$$

等效区间 $[\theta_1,\theta_2]$ 在 log 尺度上为 $[\ln 0.80,\ \ln 1.25]=[-0.2231,\ 0.2231]$。

**双单侧检验（TOST, Schuirmann 1987）**：把"等效"写成两个单侧假设

$$
H_{01}:\ \mu\le\theta_1 \quad \text{vs.}\quad H_{11}:\ \mu>\theta_1,
\qquad
H_{02}:\ \mu\ge\theta_2 \quad \text{vs.}\quad H_{12}:\ \mu<\theta_2.
$$

**两个零假设都在水平 $\alpha$ 下被拒绝**，才宣告生物等效。这等价于（Westlake / Schuirmann 的置信区间表述）：

$$
(1-2\alpha)\times100\%\ \text{置信区间}\ \subseteq\ [\theta_1,\theta_2]
\iff \text{等效成立}.
$$

> 所以 BE 常用 **90% 置信区间**（对应 $\alpha=0.05$ 的 TOST），而不是常规的 95%。

**为什么这控制了 I 类错误**：若真实比值在区间外，TOST 宣告"等效"的概率不超过 $\alpha$。

### 2.5 BE 的功效分析与样本量计算

**功效的定义**：在真实比值 $\mu$（落在 $(\theta_1,\theta_2)$ 内）处，TOST 正确宣告等效的概率：

$$
1-\beta=\Pr\!\big(T_1>t_{1-\alpha,\nu}\ \text{且}\ T_2<-t_{1-\alpha,\nu}\big).
$$

**常用的简单近似公式**（便于快速起手，但**会系统性高估功效**）。对 $2\times2$ 交叉设计、对称限 $[-\theta,\theta]$、真实比值为 $\mu$：

$$
1-\beta \;\approx\; \Phi\!\left(\frac{\sqrt{N}\,(\theta-|\mu|)}{\sigma_w\sqrt2}-t_{1-\alpha,N-2}\right),
$$

由此反解出的近似样本量为

$$
N \;\approx\; \frac{2\,\sigma_w^2\,(z_{1-\alpha}+z_{1-\beta})^2}{(\theta-|\mu|)^2}.
$$

**这个近似的缺陷**：它实际上只用了**单侧**检验的结构（相当于假设另一个单侧检验"必然通过"），忽略了 TOST 中两个检验共享同一随机量而必须**同时通过**这一事实。这个近似能把 N=12 时真实只有 64% 的功效误报成 83%，因此法规层面和高精度计算都需要**精确公式——Owen's Q**。

---

## 3. Owen's T/Q 的功效与样本量推导

### 3.1 模型与检验统计量（以 $2\times2$ 交叉设计为例）

- 总受试者数 $N$，每序列 $n=N/2$；对数尺度分析。
- 真实对数比值 $\mu$，**个体内标准差** $\sigma_w$（与个体内变异系数 CV 的关系：$\sigma_w=\sqrt{\ln(1+\text{CV}^2)}$）。
- 疗效差估计 $\hat\mu=\bar X_T-\bar X_R$，其方差（交叉设计去除了个体间变异）：

$$
\operatorname{Var}(\hat\mu)=\sigma_w^2\cdot\frac{2}{N},
\qquad
\hat\sigma_w^2:\ \frac{\nu\hat\sigma_w^2}{\sigma_w^2}\sim\chi^2_\nu,\ \ \nu=N-2.
$$

> **标准误约定说明**（易混淆点）：$2\times2$ 交叉设计的治疗差标准误是 $\mathrm{SE}=\sigma_w\sqrt{2/N}$（$N$ 为**总人数**）。有的文献写成 $\sigma_w\sqrt{2/n}$（$n$ 为每序列人数）或把 $\sigma_w$ 定义成"个体内差值的标准差"，量纲会差一个 $\sqrt2$——阅读和复现他人表格时请先确认约定，本文统一用 $\mathrm{SE}=\sigma_w\sqrt{2/N}$。

定义两个检验统计量

$$
T_1=\frac{\hat\mu-\theta_1}{\hat\sigma_w\sqrt{2/N}},\qquad
T_2=\frac{\hat\mu-\theta_2}{\hat\sigma_w\sqrt{2/N}}.
$$

**标准化**：令 $Z=\dfrac{\hat\mu-\mu}{\sigma_w\sqrt{2/N}}\sim N(0,1)$、$W=\sqrt{V/\nu}$（$V\sim\chi^2_\nu$），则

$$
T_1=\frac{Z+\delta_1}{W},\quad T_2=\frac{Z+\delta_2}{W},
\qquad
\delta_1=\frac{\mu-\theta_1}{\sigma_w\sqrt{2/N}},\quad
\delta_2=\frac{\mu-\theta_2}{\sigma_w\sqrt{2/N}}.
$$

注意两个统计量**共享同一个 $Z$ 和同一个 $W$**——这就是"相关系数为 $-1$ 的二元非中心 t"，也是 Owen (1965) 研究的情形。

### 3.2 功效作为单一积分（逐步推导）

TOST 在水平 $\alpha$ 下拒绝当且仅当

$$
T_1>t_{1-\alpha,\nu}\ \ \text{且}\ \ T_2<-t_{1-\alpha,\nu}
\iff
t_{\alpha}W-\delta_1 < Z < -t_{\alpha}W-\delta_2,
\qquad t_\alpha:=t_{1-\alpha,\nu}.
$$

这个关于 $Z$ 的区间**非空**当且仅当

$$
t_\alpha W-\delta_1 < -t_\alpha W-\delta_2
\iff
W < w^*:=\frac{\delta_1-\delta_2}{2t_\alpha}
=\frac{\theta_2-\theta_1}{2\,t_\alpha\,\sigma_w\sqrt{2/N}}.
$$

于是对 $W$ 求条件期望，得到**精确功效的单一积分**：

$$
1-\beta=\int_0^{w^*}\Big[\Phi\big(-t_\alpha w-\delta_2\big)-\Phi\big(t_\alpha w-\delta_1\big)\Big]\,f_W(w)\,dw
$$

其中 $f_W(w)=\dfrac{2(\nu/2)^{\nu/2}}{\Gamma(\nu/2)}w^{\nu-1}e^{-\nu w^2/2}$。

**对称情形**（$\mu=0$、限 $\pm\theta$，即 $\delta_1=\delta=-\delta_2$）可化简为

$$
1-\beta=\int_0^{\theta/(t_\alpha\sigma_w\sqrt{2/N})}\big[2\Phi(\delta-t_\alpha w)-1\big]\,f_W(w)\,dw,
$$

其中 $2\Phi(x)-1$ 就是"落入区间 $[-x,x]$"的概率——**两尾结构**在此体现得最清楚。

### 3.3 用 Owen's Q 表达（与单一积分等价）

利用 $\int_0^{w^*}=\int_0^{\infty}-\int_{w^*}^{\infty}$，并代入 2.3 中 $Q_\nu$ 与 $F_\nu$ 的关系，得

$$
1-\beta=F_\nu(-t_\alpha;\delta_2)-F_\nu(t_\alpha;\delta_1)-Q_\nu\!\big(-t_\alpha\sqrt\nu,\delta_2;0,R\big)+Q_\nu\!\big(t_\alpha\sqrt\nu,\delta_1;0,R\big)
$$

其中

$$
R=w^*\sqrt\nu=\frac{\sqrt\nu\,(\theta_2-\theta_1)}{2\,t_\alpha\,\sigma_w\sqrt{2/N}}.
$$

**逐项解读**：

1. $F_\nu(-t_\alpha;\delta_2)-F_\nu(t_\alpha;\delta_1)$：两个非中心 t CDF 之差——即"忽略 $W$ 的上限截断"的近似，也是 $Q_\nu(\cdot;\cdot;0,0)$ 的线性组合；
2. $-Q_\nu(-t_\alpha\sqrt\nu,\delta_2;0,R)+Q_\nu(t_\alpha\sqrt\nu,\delta_1;0,R)$：$R>0$ 处的**截断校正**。

数值验证：截断校正项在常规参数下极小，因此"两个非中心 t CDF 之差"本身就是极好的精确近似；**真正重要的**是保留两尾联合结构，而不是用单侧正态近似。

> 这也正是 Diletti et al. (1991) 与 Phillips (1990) 在 BE 文献里给出的 Owen-Q 型精确功效公式的内核。

### 3.4 样本量反推（迭代流程）

1. 用 2.5 的近似公式得到初始 $N_0$；
2. 对当前 $N$ 计算精确功效（3.2 或 3.3 式）；
3. 若功效 $<1-\beta$ 则增大 $N$，否则减小；可用二分法；
4. 取**偶数** $N$（保证两序列均衡），并常额外预留 10–20% 脱落率。

### 3.5 数值示例（工作示例）

**设定**：$2\times2$ 交叉，个体内 $\mathrm{CV}=20\%$（$\sigma_w=\sqrt{\ln 1.04}=0.19804$），真实比值 $1.0$（$\mu=0$），等效限 $0.80$–$1.25$，$\alpha=0.05$，目标功效 $80\%$。

以 $N=16$ 为例：$\nu=14$，$t_{0.95,14}=1.7613$，$\mathrm{SE}=0.19804\sqrt{2/16}=0.07002$，$\delta_1=0.22314/0.07002=3.1869$，$\delta_2=-3.1869$。代入积分得精确功效 $=0.8332\ge 0.8$；而 $N=14$ 只有 $0.7549$。因此 $N=16$。

**精确 vs 近似对比表**（同一设定）：

| $N$ | 精确功效（Owen Q） | 单侧正态近似 | "两个非中心 t CDF 之差"近似 |
|---|---|---|---|
| 12 | 0.6445 | 0.8283 | 0.6432 |
| 14 | 0.7549 | 0.8847 | 0.7548 |
| 16 | 0.8332 | 0.9230 | 0.8332 |
| 20 | 0.9249 | 0.9663 | 0.9249 |

**读表结论**：单侧正态近似系统性高估功效（$N=12$ 报 83% 而真实仅 64%），用它定样本量会导致研究**把握度不足**；"两个非中心 t CDF 之差"与精确值几乎一致，说明精确计算的核心是**保留两尾联合结构**——这正是 Owen's Q 的用途。

**其他情形的精确样本量**：

| 情形 | 目标功效 | 最小 $N$ |
|---|---|---|
| CV=20%，比值=1.0 | 80% | 16 |
| CV=20%，比值=0.95（$\mu=-0.0513$） | 80% | 19 |
| CV=20%，比值=0.95 | 90% | 25 |
| CV=30%，比值=1.0 | 80% | 32 |

> **约定提醒**：不同软件/教材由于（a）精确 vs 近似公式、（b）$\sigma_w$ 的定义、（c）80% vs 90% 功效、（d）是否含脱落补偿，给出的数字可能相差约 $\sqrt2$ 倍。本文数字对应 $\mathrm{SE}=\sigma_w\sqrt{2/N}$ 的精确 Owen-Q 计算。

---

## 4. Owen's T/Q 的数值计算方法（最新）

### 4.1 Owen's T 的计算

- **直接数值积分**：定义式在 $[0,a]$ 上做高斯求积即可，$a$ 很大时用 $T(h,\infty)=\tfrac12\Phi(-|h|)$ 截断。
- **级数/连分式（Patefield & Tandy 2000）**：把 $T(h,a)$ 表成快速收敛的级数，并用连分式加速，可稳定达到约 14–15 位有效数字。
- **经典算法**：Young & Minder (1974) 的 **AS 76**，以及 Boys (1989) 的 **AS R80** 勘误改进。
- **间接路线**：通过二元正态 CDF 算法（Genz 2004）反解 $T$。

### 4.2 非中心 t CDF 的计算

- **级数 / 不完全 beta**：Cooper (1968) **AS 5**（积分式）；Lenth (1989) **AS 243**（用正则化不完全 beta 函数的级数，是 R 的 `pt(q, df, ncp)` 的底层原型）。
- **离散混合**：Benton & Krishnamoorthy (2003) 把非中心 t 表成泊松混合的加权中心 t，对大 $\nu$、大 $|\delta|$ 更稳。
- **单一积分表示**＋高斯求积：实现简单、可任意精度。
- **渐近展开（近期）**：Gil、Segura、Temme 等近期工作给出对极端参数稳定的渐近与反演公式。

### 4.3 Owen's Q 与 TOST 功效的直接计算

- **$R=0$ 情形**：直接用非中心 t CDF。
- **$R>0$ 情形**：$Q_\nu(t,\delta;0,R)=\int_{R/\sqrt\nu}^{\infty}\Phi\!\big(\tfrac{tw}{\sqrt\nu}-\delta\big)f_W(w)\,dw$，直接高斯求积。
- **组合式**：按 3.3 式把 TOST 功效写成两个 $F_\nu$ 之差减去两个 $Q_\nu$ 校正项；实践中校正项极小。
- **高精度实现**：R 包 **`OwenQ`**（Stéphane Laurent）用 Rcpp 实现了 $Q_\nu$ 与非中心 t 的高精度版本。

### 4.4 数值稳定性要点

1. **相消问题**：$2\Phi(x)-1$ 在 $x$ 较大时应改用补误差函数 $\operatorname{erfc}$ 的尾部表达式。
2. **非中心 t 的尾部**：$|\delta|$ 很大时需用反射/互补公式。
3. **密度下溢**：$w^{\nu-1}e^{-\nu w^2/2}$ 应在对数域计算。
4. **参数边界**：$|\delta|$ 或 $\nu$ 极端时用渐近公式兜底。

---

## 5. 计算机实现（C/Fortran 公开代码）

**有公开代码**——既有 Fortran 的经典 Applied Statistics 算法，也有 C/C++ 的现代库，以及 R 生态的成熟实现：

| 库 / 算法 | 语言 | 关键函数 | 用途与来源 |
|---|---|---|---|
| AS 5 (Cooper 1968) | Fortran | 非中心 t CDF | statlib `/apstat/5` |
| AS 76 (Young & Minder 1974) | Fortran | Owen's T | statlib `/apstat/76` |
| AS R80 (Boys 1989) | Fortran | AS 76 勘误 | statlib `/apstat/R80` |
| AS 243 (Lenth 1989) | Fortran | 非中心 t CDF | statlib `/apstat/243` |
| DCDFLIB (Brown, Lovato, Russell) | C + Fortran | `cdftnc`（非中心 t）、`cdft`、`cdffnc` 等 | statlib 通用发行版 |
| Boost.Math | C++ | `owens_t`、`non_central_t`、`bivariate_normal` | boost.org（头文件库） |
| GSL | C | 仅中心 t（`gsl_cdf_tdist`） | 非中心 t 是已知空缺 |
| R `stats::pt(q, df, ncp)` | C（源自 AS 243） | 非中心 t CDF | R 基础包 |
| R `PowerTOST` (Labes & Schütz) | R/C | `power.TOST`、`sampleN.TOST`（Owen Q 精确功效） | CRAN |
| R `OwenQ` (Laurent) | R/Rcpp | `OwenQ1/Q2/Q3`、非中心 t | CRAN / GitHub |
| R `sn` (Azzalini)、`mnormt` | R | `T.Owen`、`pst`、`sadmvn`、`biv.nt.prob` | CRAN |
| SciPy `scipy.stats.nct` | C/Python | 非中心 t | SciPy |

**要点**：

- **Fortran 经典算法**（AS 系列）是"第一手"公开代码，适合直接移植或对照；DCDFLIB 同时提供 C 与 Fortran 两个版本，覆盖非中心 t。
- **现代 C++** 首选 Boost.Math：`owens_t(h, a)` 直接算 Owen's T，`non_central_t` 提供非中心 t 的 CDF/PDF/分位数，`bivariate_normal` 算二元正态。
- **BE 样本量的"开箱即用"** 首选 R 的 **`PowerTOST`**：内部正是用 Owen's Q（Phillips 1990 / Diletti 1991）做精确功效与样本量。
- **验证/高精度** 用 **`OwenQ`** 包对 $Q_\nu$ 与功效做独立复核。

---

## 6. 简要总结

1. **Owen's T/Q 是"等效/非劣"类检验功效计算的精确引擎**：它解决的问题是——当检验由**两个共享同一随机分母的非中心 t 统计量共同决定**（TOST、生物等效性）时，如何精确计算联合概率。

2. **逻辑链条**：非中心 t 的 CDF 是单一积分 $\int_0^\infty\Phi(tw-\delta)f_W(w)dw$；Owen's Q 就是这个积分的推广（下限从 $0$ 变为 $R$），且 $F_\nu(t;\delta)=Q_\nu(t\sqrt\nu,\delta;0,0)$；TOST 功效 = 两个非中心 t CDF 之差 **减** 两个 $Q_\nu$ 截断校正项，并可等价写成一个从 $0$ 到 $w^*$ 的单一积分。

3. **实践要点**：简单单侧正态近似会**系统性高估功效**（$N=12$ 时误报 83% vs 真实 64%），因此法规与高精度场景必须用精确方法；而"两个非中心 t CDF 之差"本身就已几乎精确（截断校正极小），所以精确计算的关键是**保留两尾联合结构**。

4. **数值与实现都已成熟**：Fortran 有 AS 5/76/243、DCDFLIB；C++ 有 Boost.Math；R 有 `pt(..., ncp)`、`PowerTOST`、`OwenQ`、`sn`/`mnormt`。Patefield–Tandy (2000) 与 Lenth (1989) 分别是 Owen's T 与非中心 t 计算的基准算法，近期 Gil–Segura–Temme 的渐近展开补齐了极端参数场景。

5. **一个可直接使用的精确功效公式**（$2\times2$ 交叉、$\mathrm{SE}=\sigma_w\sqrt{2/N}$、$\nu=N-2$）：

$$
1-\beta=\int_0^{w^*}\Big[\Phi\big(-t_\alpha w-\delta_2\big)-\Phi\big(t_\alpha w-\delta_1\big)\Big]\frac{2(\nu/2)^{\nu/2}}{\Gamma(\nu/2)}w^{\nu-1}e^{-\nu w^2/2}\,dw,
\qquad
w^*=\frac{\theta_2-\theta_1}{2\,t_\alpha\,\sigma_w\sqrt{2/N}}.
$$

把该积分封装成一个函数，再套用二分/迭代即可得到任意精度下的最小样本量。

---

**主要参考文献**

- Owen, D. B. (1956). Tables for computing bivariate normal probabilities. *Ann. Math. Statist.* 27, 1075–1090.
- Owen, D. B. (1965). A special case of a bivariate non-central t-distribution. *Biometrika* 52, 437–446.
- Owen, D. B. (1968). A survey of properties and applications of the non-central t-distribution. *Technometrics* 10, 445–478.
- Schuirmann, D. J. (1987). A comparison of the two one-sided tests procedure and the power approach for assessing the equivalence of average bioavailability. *J. Pharmacokinet. Biopharm.* 15, 657–680.
- Phillips, K. F. (1990). Power of the two one-sided tests procedure in bioequivalence. *J. Pharmacokinet. Biopharm.* 18, 137–144.
- Diletti, E., Hauschke, D., & Steinijans, V. W. (1991). Sample size determination for bioequivalence assessment by means of confidence intervals. *Int. J. Clin. Pharmacol. Ther. Toxicol.* 29, 1–8.
- Patefield, M., & Tandy, D. (2000). Fast and accurate computation of Owen's T function. *J. Stat. Softw.* 5(5), 1–25.
- Young, J. C., & Minder, C. E. (1974). Algorithm AS 76. *Appl. Statist.* 23, 455–457; Boys, R. J. (1989). AS R80. *Appl. Statist.* 38, 580–582.
- Lenth, R. V. (1989). Algorithm AS 243: Cumulative distribution function of the non-central t distribution. *Appl. Statist.* 38, 185–189.
- Benton, D., & Krishnamoorthy, K. (2003). Computing discrete mixtures of continuous distributions. *Comput. Statist. Data Anal.* 43, 249–267.
- Genz, A. (2004). Numerical computation of rectangular bivariate and trivariate normal and t probabilities. *Stat. Comput.* 14, 251–260.