Math in Biometrics Mathematical Analysis Leibniz's Notation Who is Leibniz 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. Leibniz's Notation System The derivative of a function based on today’s standard is given by: Leibniz's notation expresses the derivative as: where goes toward 0. Fractional Behavior Let’s review some examples where Leibniz’s notation is often utilized. The Chain Rule using Lagrange Notation states: We can translate the above Chain Rule into Leibniz's Notation as: 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. Try It on a Function Let us try . Click to Expand/Collapse Expand the function: Simplify fraction: So the derivative of is . 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 -\infty0$ 时的 $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}}. $$ 逐项解读 : $F_\nu(-t_\alpha;\delta_2)-F_\nu(t_\alpha;\delta_1)$:两个非中心 t CDF 之差——即"忽略 $W$ 的上限截断"的近似,也是 $Q_\nu(\cdot;\cdot;0,0)$ 的线性组合; $-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 样本量反推(迭代流程) 用 2.5 的近似公式得到初始 $N_0$; 对当前 $N$ 计算精确功效(3.2 或 3.3 式); 若功效 $<1-\beta$ 则增大 $N$,否则减小;可用二分法; 取 偶数 $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 数值稳定性要点 相消问题 :$2\Phi(x)-1$ 在 $x$ 较大时应改用补误差函数 $\operatorname{erfc}$ 的尾部表达式。 非中心 t 的尾部 :$|\delta|$ 很大时需用反射/互补公式。 密度下溢 :$w^{\nu-1}e^{-\nu w^2/2}$ 应在对数域计算。 参数边界 :$|\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. 简要总结 Owen's T/Q 是"等效/非劣"类检验功效计算的精确引擎 :它解决的问题是——当检验由 两个共享同一随机分母的非中心 t 统计量共同决定 (TOST、生物等效性)时,如何精确计算联合概率。 逻辑链条 :非中心 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^*$ 的单一积分。 实践要点 :简单单侧正态近似会 系统性高估功效 ($N=12$ 时误报 83% vs 真实 64%),因此法规与高精度场景必须用精确方法;而"两个非中心 t CDF 之差"本身就已几乎精确(截断校正极小),所以精确计算的关键是 保留两尾联合结构 。 数值与实现都已成熟 :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 的渐近展开补齐了极端参数场景。 一个可直接使用的精确功效公式 ($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.