搜索
关闭
搜索
高级搜索
搜索菜单

总结

考虑了多元正态分布的累积分布函数的估计。多元正态分布可以有任意正定相关矩阵和任意均值向量。所采取的方法分为两个阶段。在第一阶段，展示了如何使用递归积分方法计算非中心正交方案概率。在第二阶段，对Schläfli和Abrahamson的一些观点进行了扩展，以表明任何非中心正态概率至多可以表示为两者之间的差异(米−1)! 非中心正交方案概率。这种方法可以准确评估在统计实践中具有重要应用的许多多元正态概率。

累积分布函数,奥坦概率,正交方案概率,多面体圆锥体,递归积分

1.简介

假设随机变量x1,…,x米具有任意非奇异多元正态分布并考虑非中心正态概率

P（P）米(μ,R（右）)=公共关系x我⩾0,1⩽我⩽米=∫0∞…∫0∞ϕ米(x;μ,R（右）)d日x1…d日x米
1

哪里x=(x1,…,x米)′∼N个米(μ,R（右）)带平均向量μ=(μ1,…,μ米)′和正定相关矩阵R（右）={ρij公司}具有密度函数φ米(x;μ,R（右）). 注意相关性矩阵R（右）不更改，如果x替换为−x=(−x1,…,−x米)'多元正态分布函数可以表示为

F类米(c)=公共关系x我⩽c我,1⩽我⩽米=公共关系-x我⩾-c我,1⩽我⩽米=P（P）米(-μ+c,R（右）)

对于任何向量c=(c1,…,c米)′.

由于多元正态分布在统计应用中起着基础性作用，关于积分表达式（1）的计算已有大量文献。除非尺寸米很小，任何通过数值积分技术直接计算表达式（1）的方法都是不可管理的。替代的近似方法包括无限级数展开、数值积分中的降维和蒙特卡罗模拟。有关此问题的一般讨论，请参见，例如，约翰逊和科茨（1972），第35章，斯图亚特和奥德（1994）第15章，以及童（1990），第8章。

如果μ=0=（0，…，0）′，则P（P）米(0,R（右）)是中心正态概率，此外，对于三对角相关矩阵R（右）令人满意的ρij公司=0用于|我−j个| > 1,P（P）米(μ,R（右）)称为非中心正交概率P（P）米(0,R（右）)是中心正交概率。正态概率和正态概率之间有一些有趣的关系。Schläfli（1858）证明了中心正态概率可以写成米! 中心正交方案概率。亚伯拉罕森（1964）用米=4四元中心正态概率可以表示为不超过六个（而不是4！=24）中心正态分布概率的线性组合。

本文推广了这些结果，并说明了如何将任何非中心正态概率表示为至多(米−1)! 非中心正交方案概率。此外，还提出了一种快速准确的递归积分算法来计算非中心正交方案的概率。总之，这些结果提供了一种准确评估表达式（1）以及任何多元分布函数的方法F类米(c).

首先，在第节2给出了计算非中心正交方案概率的递推积分方法。然后，如第节所述三，所需的非中心正态概率用非中心正交方案概率表示。论文以第节结束4其中给出了该方法的一些数值示例。

2.非中心正交方案概率的评估

在本节中，关注的是米×米三对角相关矩阵R（右）={ρij公司}为此ρij公司=0（对于）|我−j个| > 1.三輪等. (2000)给出了计算中心正交方案概率的过程P（P）米(0,R（右）). 然而，在本节中，提出了一种改进的程序，用于评估任何非中心正交方案的概率P（P）米(μ,R（右）).

2.1. 递归积分公式

正定三对角矩阵R（右）可以分解为R（右）=B类B类'，其中的唯一非零元素B类={b条ij公司}是对角元素b条ii（ii）,1≤我≤米和对角线下方的元素b条我,我−1,2≤我≤米。如果我们定义

D类我=1ρ210⋯0ρ211ρ32⋯0⋮⋱⋱⋱⋮0⋯ρ我-1,我-21ρ我,我-10⋯0ρ我,我-111⩽我⩽米,

然后D类我> 0,1≤我≤米，自R（右）是正定的。此外，这些行列式由方程式关联

D类我=D类我-1-ρ我,我-12D类我-2,2⩽我⩽米,

具有D类0=D类1=1，因此可以按顺序计算。此外B类然后可以通过将R（右）=B类B类'作为

b条11=1=D类1/D类0,b条我,我-1=ρ我,我-1/b条我-1,我-1=ρ我,我-1D类我-2/D类我-1,b条我我=1-b条我,我-12=D类我/D类我-1,2⩽我⩽米.

下一步，如果x=(x1,…,x米)′∼N个米(μ,R（右）)，然后根据分解R（右）=B类B类“那个z（z）=(z（z）1,…,z（z）米)′=B类−1(x−μ)∼N个米(0,我米)和

x1=z（z）1+μ1,x我=b条我,我-1z（z）我-1+b条我我z（z）我+μ我,2⩽我⩽米.

然后

P（P）米(μ,R（右）)=公共关系x我⩾0,1⩽我⩽米=公共关系z（z）1+μ1⩾0,b条我,我-1z（z）我-1+b条我我z（z）我+μ我⩾0,2⩽我⩽米

现在定义

（f）米-1(z（z）)=∫-μ米-b条米,米-1z（z）/b条米米∞ϕ(t吨)d日t吨,（f）我-1(z（z）)=∫-μ我-b条我,我-1z（z）/b条我我∞（f）我(t吨)ϕ(t吨)d日t吨,2⩽我⩽米-1,
2

哪里φ(t吨)是标准正态概率密度函数，因此所需概率由下式给出

P（P）米(μ,R（右）)=∫-μ1∞（f）1(z（z）)ϕ(z（z）)d日z（z）.
三

这些积分表达式可以使用如下递归计算方法进行计算。首先，功能（f）米−1(z（z）)在精细的点网格上进行计算，并且值存储为数组。然后使用方程式（2）的值（f）米−2(z（z）)可以在类似的点网格上计算，也可以存储为数组。此时可以删除原始数组。继续以这种方式，功能（f）我(z（z）)通过递归计算，最终通过计算方程（3）。计算强度在米因此，对于任何实际值米通过合并足够数量的网格点。

这种递归计算方法已用于Hayter和Liu（1996）和三輪等. (2000)也应该指出莫兰（1985）当ρ我,我−1都等于某些值x1,…,x米依据米+1个独立的正态随机变量。

2.2. 与精确值的比较

什么时候？ρ我,我−1=ρ,2≤我≤米、和ρij公司=0,|我−j个|>1，精确值P（P）米(0,R（右）)已知的一些特殊值ρ例如，如果ρ=-12然后P（P）米(0,R（右）)=1/(米+1)!. 此外，当ρ=12莫兰（1983）证明了

P（P）2n个(0,R（右）)=(-1)n个×22n个+222n个+2-1B类2n个+2(2n个+2)!,P（P）2n个+1(0,R（右）)=(-1)n个+1E类2n个+2(2n个+2)!,

哪里B类n个和E类n个是伯努利数和欧拉数。

这些结果允许检查建议的递归积分过程的准确性。在递归积分（2）的实现中，网格点在[-8,8]内分配，其中网格间隔的长度与φ(z（z）)−1/4和三次多项式近似应用于函数（f）我(z（z）). 网格点的选择和三次多项式近似在中进行了详细讨论三輪等. (2000)。的结果米=5和米=10，见表1哪里G公司是使用的网格点数量。使用G公司=128个网格点，数值结果与已知精确值一致，小数点后八位，但情况除外米=10和ρ=-12⁠。然而，使用G公司=案例的512个网格点米=10和ρ=-12⁠，其中概率值实际上非常小。在所有情况下，计算时间都可以忽略不计。

表1
新标签中打开

正交方案概率的数值计算†

G公司m和以下值的概率ρ:
米=5米=10
ρ=12ρ=-12ρ=12ρ=-12
320.08472210.0013888850.008863210.2502×10−7
640.0847222190.00138888880.00886323450.25051×10−7
1280.0847222220.00138888890.00886323550.250520×10−7
2560.25052105×10−7
5120.25052108×10−7
完全正确0.0847222220.00138888890.00886323550.25052108×10−7
G公司m和以下值的概率ρ:
米=5米=10
ρ=12ρ=-12ρ=12ρ=-12
320.08472210.0013888850.008863210.2502×10−7
640.0847222190.00138888880.00886323450.25051×10−7
1280.0847222220.00138888890.00886323550.250520×10−7
2560.25052105×10−7
5120.25052108×10−7
完全正确0.0847222220.00138888890.00886323550.25052108×10−7
†

G公司是网格点的数量。

表1
新标签中打开

正交方案概率的数值计算†

G公司m和以下值的概率ρ:
米=5米=10
ρ=12ρ=-12ρ=12ρ=-12
320.08472210.0013888850.008863210.2502×10−7
640.0847222190.00138888880.00886323450.25051×10−7
1280.0847222220.00138888890.00886323550.250520×10−7
2560.25052105×10−7
5120.25052108×10−7
完全正确0.0847222220.00138888890.00886323550.25052108×10−7
G公司m和以下值的概率ρ:
米=5米=10
ρ=12ρ=-12ρ=12ρ=-12
320.08472210.0013888850.008863210.2502×10−7
640.0847222190.00138888880.00886323450.25051×10−7
1280.0847222220.00138888890.00886323550.250520×10−7
2560.25052105×10−7
5120.25052108×10−7
完全正确0.0847222220.00138888890.00886323550.25052108×10−7
†

G公司是网格点的数量。

应该记住，对于非中心正交方案概率，递归积分过程不会增加计算难度P（P）米(μ,R（右）)，尽管在这种情况下没有闭合形式的表达式可用于概率。第节讨论了评估多元正态概率的整体方法的准确性4.

3.将正态概率表示为正交方案概率的差异

在本节中，我们将展示非中心正态概率P（P）米(μ,R（右）)最多可以表示为(米−1)! 非中心正交方案概率。这扩展了Schläfli（1858）并基于亚伯拉罕森（1964）具有米=4.

为了继续，有以下关于对称矩阵形式的定义是有用的Σ.

定义1。对于0≤第页≤米−3安米×米对称矩阵Σ={σij公司}据说有“正交序”第页'如果σij公司1≤时=0我≤第页,我+2≤j个≤米和σ第页+1,j个≠0对于某些j个令人满意的第页+3≤j个≤米.三对角矩阵定义为具有正交方案阶米−2.

根据这个定义，任何具有σ1j个≠0对于某些j个满足3≤j个≤米具有正交方案阶0。然而，如果正态随机变量x1,…,x米具有协方差矩阵Σ，则与这些随机变量的某些排列相对应的协方差矩阵可能具有不同于Σ.

本节的主要结果是展示如何根据相关矩阵定义非中心正态概率R（右）正交方案阶第页最多可以用米−第页−1非中心正态概率定义为相关矩阵，其正交方案阶严格大于第页.

正定相关矩阵R（右）可以写在表格中R（右）=A类′A类其中矩阵A类由线性无关的列向量组成一1,…,一米，所以一我′一我=1和一我′一j个=ρij公司对于我≠j个.让多面体圆锥体问，顶点位于原点，由定义

问=z（z）:一我'z（z）⩾0,1⩽我⩽米.

那么如果x∼N个米(μ,R（右）)因此z（z）=A类′−1x∼N个米(A类′−1μ,我米)，因此

P（P）米(μ,R（右）)=公共关系x我⩾0,1⩽我⩽米=公共关系一我'z（z）⩾0,1⩽我⩽米=公共关系{z（z）∈问}.

多面体圆锥体问以为界米超平面

H（H）我=z（z）:一我'z（z）=0,1⩽我⩽米,

哪里一我是超平面的法向量H（H）我此外，了解圆锥体的边向量也很有用。假设五1,…,五米是的列向量V（V）=A类′−1以便一我′五我=1和一我′五j个=0（对于）我≠j个.矢量五k个同时属于超平面H（H）我,1≤我≤米,我≠k个，因此可以将其视为多面体圆锥体的边向量问.多面体圆锥体问然后可以用线性无关的边向量表示五1,…,五米作为

问=z（z）:z（z）=λ1五1+…+λ米五米,λ我⩾0.

如果我们考虑c我五我而不是五我对于任何正常数c我或者如果向量五1,…,五米以不同的顺序重新排列。

现在假设相关矩阵R（右）={ρij公司}具有正交序列第页0≤第页≤米−3.让一个新的边向量第页由定义

第页=γ第页+2五第页+2+…+γ米五米,

其中系数γ秒,第页+2≤秒≤米，由给出

γ秒=-ρ第页+1,秒ρ第页+1,秒如果有的话ρ第页+1,我⩽0的第页+2⩽我⩽米, 否则。 
4

所以至少有一个γ秒是绝对肯定的。接下来，针对每个γ秒≠0,第页+2≤秒≤米，考虑一个圆锥体问(秒)由线性无关的边缘向量定义

五1,…,五秒-1,第页,五秒+1,…,五米.
5

然后我们可以建立以下引理。证据见附录A.

引理1

  • （a）

    问∪(∪γ秒<0问(秒))=∪γ秒>0问(秒).

  • （b）

    十字路口问(秒)∩问(t吨),秒≠t吨，包含在维度的子空间中米−1，如果γ秒和γt吨有相同的标志。

  • （c）

    十字路口问∩问(秒)对于γ秒<0包含在维度的子空间中米−1.

由于维度的每个子空间米−1的概率为0，它直接从引理1得出

公共关系{z（z）∈问}+∑γ秒<0公共关系z（z）∈问(秒)=∑γ秒>0公共关系z（z）∈问(秒).

此外，通过按顺序重新排列边缘向量（5）

五1,…,五第页+1,第页,五第页+2,…,五秒-1,五秒+1,…,五米,

相应的法向量一我(秒),1⩽我⩽米⁠，被视为

一我(秒)=一我,1⩽我⩽第页+1,一第页+2(秒)=sgn公司γ秒一秒,一我(秒)=c我-1(秒)一我-1-ρ第页+1,我-1ρ第页+1,秒一秒,第页+三⩽我⩽秒,一我(秒)=c我(秒)一我-ρ第页+1,我ρ第页+1,秒一秒,秒+1⩽我⩽米,

具有

c我(秒)=1-2ρ第页+1,我ρ我秒ρ第页+1,秒+ρ第页+1,我2ρ第页+1,秒2-1/2.

请注意一我(秒)/一j个(秒)=0对于1≤我≤第页+1,我+2≤j个≤米然后定义R（右）(秒)=ρ我j个(秒)和μ(秒)=μ1(秒),…,μ米(秒)'通过

ρ我j个(秒)=一我(秒)'一j个(秒),μ我(秒)=一我(秒)'A类'-1μ
6

以便

公共关系z（z）∈问(秒)=公共关系一我(秒)'z（z）⩾0,1⩽我⩽米=P（P）米μ(秒),R（右）(秒),

这一部分的结果总结在以下定理中。

定理1

假设正定相关矩阵R（右）={ρij公司}具有正交序列第页0≤第页≤米−3.然后，对于任何矢量μ,

P（P）米(μ,R（右）)=∑γ秒>0P（P）米μ(秒),R（右）(秒)-∑γ秒<0P（P）米μ(秒),R（右）(秒),
7

哪里γ秒,第页+2≤秒≤米，由表达式（4）和R（右）(秒)和μ(秒)通过表达式（6）。

定理1的要点是相关矩阵R（右）(秒)每个矩阵都有一个严格大于相关矩阵的正交级数的正交级数R（右）具体来说，如果R（右）具有正交序列米−3然后是矩阵R（右）(秒)都是三对角矩阵。因为最多有米−第页−1个非零元素γ秒,第页+2≤秒≤米根据定理1的重复应用，可以得出任何非中心正态概率P（P）米(μ,R（右）)最终最多可以分解为(米−1)! 非中心正交方案概率。还应指出R（右）(秒)和μ(秒)定理1中给出了一个执行这种分割的实用算法。

第节中描述的递归积分方法2允许快速准确地评估任意值的正交方案概率米.与最多涉及的解剖同时进行(米−1)! 本节中给出的正态项，因此提出的程序允许准确评估非中心正态概率P（P）米(μ,R（右）)对于任何相关矩阵R（右）.

4.数值结果

自从第节中给出的解剖三表示非中心正态概率P（P）米(μ,R（右）)由于非中心正交方案概率之间的差异是精确的，因此该程序计算多元正态概率的数值精度仅取决于正交方案概率计算的数值精度。为了研究程序的整体数值精度，本节介绍了一些计算方法，用于通过其他方法评估正值概率的准确值。这些计算是在一台配备1.5 GHz Intel®Pentium®4处理器的个人计算机上用C语言实现的。

4.1. 中心正值概率

首先考虑相关矩阵R（右）具有ρ我j个=12,我≠j个⁠在这种情况下，众所周知P（P）米(0,R（右）)=1/(米+1). 此外，考虑协方差矩阵Σ其逆函数Σ−1={σij公司}具有个元素

σ我我=1,σ我,我+1=σ我+1,我=-12,σ我j个=0,|我-j个|>1;

在这种情况下阿尼斯和劳埃德（1953）证明了相应的中心正态概率P（P）米(0,R（右）)也等于1/(米+1). 对应的相关矩阵Σ所有条目都为非零，但当节中的解剖三所需正交方案概率的数量大大小于(米−1)!. 例如，对于米=9与（9−1）相比=40 320. 为了研究累加正交方案概率对总体精度的影响，这两个中心正态概率是通过针对这种情况提出的程序计算的米=9，其中已知准确概率为0.1。结果如所示表2可以看出，在这两种情况下，获得的值与准确值非常接近。

表2
新标签中打开

正值概率的数值计算†

G公司以下值的概率ρ或σ:
ρ我j个=12σ我我=1,σ我,我+1=-12
320.0999751100 (12.61)0.1000142803 (0.08)
640.0999990004 (23.01)0.1000006159 (0.16)
1280.0999999489 (44.02)0.1000000321 (0.30)
2560.0999999968 (87.49)0.1000000020 (0.59)
5120.0999999998 (172.12)0.1000000001 (1.18)
完全正确0.10000000000.1000000000
G公司以下值的概率ρ或σ:
ρ我j个=12σ我我=1,σ我,我+1=-12
320.0999751100 (12.61)0.1000142803 (0.08)
640.0999990004 (23.01)0.1000006159 (0.16)
1280.0999999489 (44.02)0.1000000321 (0.30)
2560.0999999968 (87.49)0.1000000020 (0.59)
5120.0999999998 (172.12)0.1000000001 (1.18)
完全正确0.10000000000.1000000000
†

G公司是网格点的数量。括号中给出了以秒为单位的计算时间。

表2
新标签中打开

正值概率的数值计算†

G公司以下值的概率ρ或σ:
ρ我j个=12σ我我=1,σ我,我+1=-12
320.0999751100 (12.61)0.1000142803 (0.08)
640.0999990004 (23.01)0.1000006159 (0.16)
1280.0999999489 (44.02)0.1000000321 (0.30)
2560.0999999968 (87.49)0.1000000020 (0.59)
5120.0999999998 (172.12)0.1000000001 (1.18)
完全正确0.10000000000.1000000000
G公司以下值的概率ρ或σ:
ρ我j个=12σ我我=1,σ我,我+1=-12
320.0999751100 (12.61)0.1000142803 (0.08)
640.0999990004 (23.01)0.1000006159 (0.16)
1280.0999999489 (44.02)0.1000000321 (0.30)
2560.0999999968 (87.49)0.1000000020 (0.59)
5120.0999999998 (172.12)0.1000000001 (1.18)
完全正确0.10000000000.1000000000
†

G公司是网格点的数量。括号中给出了以秒为单位的计算时间。

4.2. 非中心正值概率

非中心正值概率的附加计算P（P）米(μ,R（右）)对与ρij公司=ρ,我≠j个、和μ我=μ,1≤我≤米.剖宫产手术时三可以看出(米−1)! 需要正交方案概率，但它们都是相同的。这些正态概率也可以通过涉及标准正态累积分布函数的项的一维积分来计算。在童（1990）对于米=2,3,…,10,μ=−2.0、−1.0、…、4.0和ρ=0.1、0.2、…、0.9（以及许多其他情况）。每一个正态概率都是通过建议的程序计算出来的G公司=128个网格点，并与Tong的表格进行了比较。几乎所有的9×7×9=567条目都完全一致，其中有11条条目在小数点后五位相差1个单位，有两条条目在第五位相差2个单位。然而，最后两个条目与表中的值一致古普塔（1963）(表3).

表3
新标签中打开

正值概率的比较

参数以下来源的概率：
拟议程序童（1990）古普塔（1963）
米=9,ρ=0.9,μ=0.00.313800.313820.31380
米=10,ρ=0.9,μ=0.00.307470.307490.30747
参数以下来源的可能性：
拟议程序童（1990）古普塔（1963）
米=9,ρ=0.9,μ=0.00.313800.313820.31380
米=10,ρ=0.9,μ=0.00.307470.307490.30747
表3
新标签中打开

正值概率的比较

参数以下来源的可能性：
拟议程序童（1990）古普塔（1963）
米=9,ρ=0.9,μ=0.00.313800.313820.31380
米=10,ρ=0.9,μ=0.00.307470.307490.30747
参数以下来源的可能性：
拟议程序童（1990）古普塔（1963）
米=9,ρ=0.9,μ=0.00.313800.313820.31380
米=10,ρ=0.9,μ=0.00.307470.307490.30747

总的来说，这些计算表明，本文提出的方法对于有中心和无中心正态概率的数值精度都是令人满意的。

4.3. 关于计算时间的一点注记

计算正交方案概率的计算时间与米×G公司哪里G公司是第节中描述的递归方法中的网格点数量2因此，如果有适当的概率P（P）米(μ,R（右）)被分解成(米−1)! 正交方案概率，总计算时间与米!×G公司.表2证实了计算时间在G公司。也可以从中观察到表1和2通过加倍网格点的数量，可以获得至少多一个小数位的精度。

研究了不同值的计算时间米在等相关情况下ρ我j个=12其中(米−1)! 需要正交方案概率(表4). 可以看出，计算时间以米! 因此，计算将变得难以进行米然而，如果相关矩阵具有特殊结构，那么三产生了一些小于(米−1)!, 那么我们的程序可能仍然适用于米.

表4
新标签中打开

已知精确概率的计算时间比较†

米我们的程序结果和以下G值：Genz程序的结果和以下值ɛ:
G公司=16G公司=128ɛ=0.001ɛ=0.0001
时间（s）绝对误差时间（s）绝对误差时间（s）绝对误差时间（s）绝对误差
502.1×10−60.011.4×10−90.682.1×10−459.673.5×10−5
60.013.5×10−60.082.0×10−100.742.7×10−497.862.2×10−5
70.079.2×10−60.585.8×10−91.143.5×10−4100.222.7×10−5
80.585.7×10−54.762.1×10−81.102.5×10−497.703.3×10−5
95.351.6×10−444.025.1×10−81.552.8×10−4137.623.6×10−5
1053.993.5×10−4450.461.0×10−71.392.1×10−4124.862.3×10−5
201.852.8×10−4198.903.3×10−5
米我们的程序结果和以下G值：Genz程序的结果和以下值ɛ:
G公司=16G公司=128ɛ=0.001ɛ=0.0001
时间（s）绝对误差时间（s）绝对误差时间（s）绝对误差时间（s）绝对误差
502.1×10−60.011.4×10−90.682.1×10−459.673.5×10−5
60.013.5×10−60.082.0×10−100.742.7×10−497.862.2×10−5
70.079.2×10−60.585.8×10−91.143.5×10−4100.222.7×10−5
80.585.7×10−54.762.1×10−81.102.5×10−497.703.3×10−5
95.351.6×10−444.025.1×10−81.552.8×10−4137.623.6×10−5
1053.993.5×10−4450.461.0×10−71.392.1×10−4124.862.3×10−5
201.852.8×10−4198.903.3×10−5
†

等相关相关矩阵：ρ我j个≡ρ=12,μ我=0.0⁠.G公司是网格点的数量。的价值ɛGenz的程序使用它来停止重复(Genz，1992年).

表4
新标签中打开

已知精确概率的计算时间比较†

米我们的程序结果和以下G值：Genz程序的结果和以下值ɛ:
G公司=16G公司=128ɛ=0.001ɛ=0.0001
时间（s）绝对误差时间（s）绝对误差时间（s）绝对误差时间（s）绝对误差
502.1×10−60.011.4×10−90.682.1×10−459.673.5×10−5
60.013.5×10−60.082.0×10−100.742.7×10−497.862.2×10−5
70.079.2×10−60.585.8×10−91.143.5×10−4100.222.7×10−5
80.585.7×10−54.762.1×10−81.102.5×10−497.703.3×10−5
95.351.6×10−444.025.1×10−81.552.8×10−4137.623.6×10−5
1053.993.5×10−4450.461.0×10−71.392.1×10−4124.862.3×10−5
201.852.8×10−4198.903.3×10−5
米我们的程序结果和以下G值：Genz程序的结果和以下值ɛ:
G公司=16G公司=128ɛ=0.001ɛ=0.0001
时间（s）绝对误差时间（s）绝对误差时间（s）绝对误差时间（s）绝对误差
502.1×10−60.011.4×10−90.682.1×10−459.673.5×10−5
60.013.5×10−60.082.0×10−100.742.7×10−497.862.2×10−5
70.079.2×10−60.585.8×10−91.143.5×10−4100.222.7×10−5
80.585.7×10−54.762.1×10−81.102.5×10−497.703.3×10−5
95.351.6×10−444.025.1×10−81.552.8×10−4137.623.6×10−5
1053.993.5×10−4450.461.0×10−71.392.1×10−4124.862.3×10−5
201.852.8×10−4198.903.3×10−5
†

等相关相关矩阵：ρ我j个≡ρ=12,μ我=0.0⁠.G公司是网格点的数量。的价值ɛGenz的程序使用它来停止重复(Genz，1992年).

另一种可能的正态概率计算方法是蒙特卡罗方法。Genz（1992）提出了一种结合积分变量变换的蒙特卡罗方法。他的方法是针对同一组正值概率进行评估的(表4). 结果表明，Genz方法的计算时间不超过几秒米大到20，只要绝对误差小于10−3然而，蒙特卡罗方法的一个主要困难是计算时间增加了大约100倍，以实现精度的小数点后一位的增加。在精确概率未知的情况下，还比较了随机生成的相关矩阵和平均向量的计算时间(表5). 对于中等值米（例如。米≤7). 因此，我们可以期望通过增加网格点的数量可以获得足够的精度G公司，因为计算时间在G公司.

表5
新标签中打开

未知准确概率的计算时间比较†

米我们程序的结果‡,ɛ=0.001Genz手术结果,ɛ=0.001
计算时间（s）绝对误差§(×10−4)计算时间（s）绝对误差§(×10−4)
平均§§标准偏差§平均§§标准偏差§§平均§§标准偏差§§平均§§标准偏差§§
50.040.042.47.15.623.802.81.7
60.580.262.43.78.155.503.12.1
74.541.872.22.89.155.202.92.3
833.389.8859.812.707.983.42.9
米我们程序的结果‡,ɛ=0.001Genz程序的结果,ɛ=0.001
计算时间（s）绝对误差§(×10−4)计算时间（s）绝对误差§(×10−4)
平均§§标准偏差§平均§§标准偏差§§平均§§标准偏差§§平均§§标准偏差§§
50.040.042.47.15.623.802.81.7
60.580.262.43.78.155.503.12.1
74.541.872.22.89.155.202.92.3
833.389.8859.812.707.983.42.9
†

随机相关矩阵和随机平均向量。

‡

我们的程序从G公司=16个网格点。通过加倍网格点的数量重复计算G公司直到两个连续计算概率的绝对差值小于ε/2.

§

绝对误差是根据我们的程序计算的正态概率得出的G公司=1024个网格点。

§§

通过20组随机生成的相关矩阵和平均向量计算平均和标准偏差SD。

表5
新标签中打开

未知准确概率的计算时间比较†

米我们程序的结果‡,ɛ=0.001Genz程序的结果,ɛ=0.001
计算时间（s）绝对误差§(×10−4)计算时间（s）绝对误差§(×10−4)
平均§§标准偏差§平均§§标准偏差§§平均§§标准偏差§§平均§§标准偏差§§
50.040.042.47.15.623.802.81.7
60.580.262.43.78.155.503.12.1
74.541.872.22.89.155.202.92.3
833.389.8859.812.707.983.42.9
米我们程序的结果‡,ɛ=0.001Genz程序的结果,ɛ=0.001
计算时间（s）绝对误差§(×10−4)计算时间（s）绝对误差§(×10−4)
平均§§标准偏差§平均§§标准偏差§§平均§§标准偏差§§平均§§标准偏差§§
50.040.042.47.15.623.802.81.7
60.580.262.43.78.155.503.12.1
74.541.872.22.89.155.202.92.3
833.389.8859.812.707.983.42.9
†

随机相关矩阵和随机平均向量。

‡

我们的程序从G公司=16个网格点。通过加倍网格点的数量重复计算G公司直到两个连续计算概率的绝对差值小于ε/2.

§

绝对误差是根据我们的程序计算的orthant概率获得的G公司=1024个网格点。

§§

通过20组随机生成的相关矩阵和平均向量计算平均和标准偏差SD。

5.结束语

我们已经展示了非中心正态概率P（P）米(μ,R（右）)=优先级{x我≥0,1≤我≤米}多元正态随机变量x∼N个米(μ,R（右）)可以对任何相关矩阵进行准确评估R（右）和任何平均向量μ这使我们能够计算多元正态分布函数F类米(c)=优先级{x我≤c我,1≤我≤米}. 此外，已知双边概率

F类米(c,d日)=公共关系c我⩽x我⩽d日我,1⩽我⩽米

可以用2表示米非中心正值概率。

提出的程序基于递归积分算法来评估非中心正交方案概率，并且我们已经证明，利用相关矩阵的三对角性质，计算强度在维数上仅为线性米同样有趣的是，双边概率F类米(c,d日)对于三对角矩阵R（右）也可以通过类似于第节所述的递归算法直接进行计算2.

最后，用C语言编写了实现所建议程序的计算机程序，可从第一作者处获得。

致谢

作者感谢联合主编、副主编和裁判的宝贵意见和建议。

工具书类

亚伯拉罕森
,
I.G.公司。
(
1964
)
四元正态分布的Ortant概率
.
安。数学。统计师。
,
35
,
1685
–
1703
.

谷歌学者

交叉参考
搜索ADS

 

阿尼斯
,
答：A。
和
劳埃德
,
E.H.公司。
(
1953
)
关于有限个独立正态变量的部分和的范围
.
生物特征
,
40
,
35
–
42
.

谷歌学者

交叉参考
搜索ADS

 

根茨
,
A。
(
1992
)
多元正态概率的数值计算
.
J.计算。图表。统计师。
,
1
,
141
–
149
.

谷歌学者

OpenURL占位符文本

 

古普塔
,
S.S.公司。
(
1963
)
多元正态和多元的概率积分t吨
.
安。数学。统计师。
,
34
,
792
–
828
.

谷歌学者

交叉参考
搜索ADS

 

海特
,
A.J.公司。
和
线路接口单元
,
西。
(
1996
)
关于&pr；计算的注释；{X（X）1<X（X）2< ⋯ <X（X）k个}
.
美国统计局
,
50
,
365
.

谷歌学者

OpenURL占位符文本

 

约翰逊
,
N.编号。
和
克茨
,
美国。
(
1972
)
统计学中的分布：连续多元分布
.
纽约
:
威利
.

谷歌学者

OpenURL占位符文本

 

三輪
,
T。
,
海特
,
A.J.公司。
和
线路接口单元
,
西。
(
2000
)
不等方差正态随机变量水平概率的计算及其在非平衡单向模型Bartholomew检验中的应用
.
计算。统计师。数据分析。
,
34
,
17
–
32
.

谷歌学者

交叉参考
搜索ADS

 

莫兰
,
P.A.P.公司。
(
1983
)
多元正态分布的一种新的展开式
.
澳大利亚。J.统计。
,
25
,
339
–
344
.

谷歌学者

交叉参考
搜索ADS

 

莫兰
,
P.A.P.公司。
(
1985
)
多元正态概率的计算——另一种特殊情况
.
澳大利亚。J.统计。
,
27
,
60
–
67
.

谷歌学者

交叉参考
搜索ADS

 

施拉弗利
,
L。
(
1858
)
关于重积分n个dxdy公司…第纳尔，其限制为第页1=一1x+b条1年+…+小时1z（z）> 0,第页2> 0,…,第页n个>0和x2+年2+…+z（z）2< 1
.
Q.J.纯应用。数学。
,
2
,
269
–
301
.

谷歌学者

OpenURL占位符文本

 

斯图亚特
,
A。
和
订单
,
J.K。
(
1994
)
肯德尔的高级统计理论
，卷。
1
第5版。
伦敦
:
阿诺德
.

谷歌学者

OpenURL占位符文本

 

用钳子钳起
,
Y.L.公司。
(
1990
)
多元正态分布
.
纽约
:
施普林格
.

谷歌学者

交叉参考
搜索ADS

 

附录A：引理1的证明

A.1、。零件证明（a）

考虑任何λ我≥0,1≤我≤米，假设minγ秒>0(λ秒/γ秒)=λt吨/γt吨具有γt吨> 0. 然后z（z）=λ1五1+…+λ米五米∈问可以写为

z（z）=λ1五1+…+λt吨五t吨+…+λ米五米=λ1五1+…+λt吨γt吨第页-∑我=第页+2,我≠t吨米γ我五我+…+λ米五米

式中，系数第页是λt吨/γt吨≥0，系数五我,1≤我≤第页+1，是λ我≥0和系数五我,第页+2≤我≤米,我≠t吨，是λ我−λt吨γ我/γt吨≥0。因此z（z）∈问(t吨)具有γt吨>0，这意味着问⊆∪γ秒>0问(秒)此外，

z（z）=λ1五1+…+λ秒第页+…+λ米五米∈问(秒)

对于γ秒<0可以写为

z（z）=λ1五1+…+λ秒第页+…+λt吨γt吨第页-∑我=第页+2,我≠t吨米γ我五我+…+λ米五米

式中，系数第页是λ秒+λt吨/γt吨≥0，系数五我,1≤我≤第页+1，是λ我≥0，系数五秒是−λt吨γ秒/γt吨≥0和系数五我,第页+2≤我≤米,我≠秒,我≠t吨，是λ我−λt吨γ我/γt吨≥0。因此，z（z）∈问(t吨)具有γt吨>0，这意味着γ秒<0问(秒)⊆∪γ秒>0问(秒).

此外，

z（z）=λ1五1+…+λ秒第页+…+λ米五米∈问(秒)

具有λ我≥0,1≤我≤米，用于γ秒>0可以写为

z（z）=λ1五1+…+λ秒∑我=第页+2米γ我五我+…+λ米五米.

在这个表达式中五我,1≤我≤第页+1，是λ我≥0，系数五秒是λ秒γ秒≥0五我,第页+2≤我≤米,我≠秒，是λ我+λ秒γ我.让W公司是一组索引我,第页+2≤我≤米,我≠秒，该系数严格为负。如果W公司=∅，则z（z）∈问但是，如果W公司≠∅注意γ我<0用于我∈W公司假设max我∈W公司(λ我/γ我)=λw个/γw个具有γw个< 0. 然后z（z）可以表示为

z（z）=λ1五1+…+λ秒第页+…+λw个五w个+…+λ米五米=λ1五1+…+λ秒第页+…+λw个γw个第页-∑我=第页+2,我≠w个米γ我五我+…+λ米五米,
8

式中，系数第页是λ秒+λw个/γw个=(λw个+λ秒γw个)/γw个≥0，系数五我,1≤我≤第页+1，是λ我≥0，系数五秒是−λw个γ秒/γw个≥0且系数为五我,第页+2≤我≤米,我≠秒,我≠w个，是λ我−λw个γ我/γw个.如果γ我≥0则清晰λ我−λw个γ我/γw个≥0，如果我∈W公司那么根据以下定义，系数也同样为非负w个最后，如果γ我<0和我∉W公司然后，通过将系数写为

λ我-λw个γ我γw个=λ我+λ秒γ我-γ我γw个λw个+λ秒γw个⩾0.

因此，表达式（8）中的所有系数z（z）非负，这意味着如果W公司≠∅则z（z）∈问(w个)具有γw个<0.因此，可以证明γ秒>0问(秒)⊆问∪(∪γ秒<0问(秒))而且已经确定

问∪∪γ秒>0问(秒)=∪γ秒>0问(秒).

A.2、。零件证明（b）

假设γ秒和γt吨,秒≠t吨，有相同的符号。如果z（z）∈问(秒)∩问(t吨)它可以表示为

z（z）=λ1五1+…+λ秒第页+…+λ米五米=λ1五1+…+λ秒γ第页+2五第页+2+…+γ米五米+…+λ米五米=ξ1五1+…+ξt吨第页+…+ξ米五米=ξ1五1+…+ξt吨γ第页+2五第页+2+…+γ米五米+…+ξ米五米,

哪里λ我≥0且ξ我1≤时≥0我≤米.等式系数五秒和五t吨给予

λ秒γ秒=ξ秒+ξt吨γ秒,λt吨+λ秒γt吨=ξt吨γt吨,

因此λt吨γ秒=−ξ秒γt吨。自γ秒和γt吨非零且符号相同，这要求λt吨=ξ秒=0。因此，问(秒)∩问(t吨)包含在维的子空间中米−1.

答3。零件证明（c）

如果z（z）∈问∩问(秒)对于γ秒<0可以表示为

z（z）=λ1五1+…+λ秒五秒+…+λ米五米=ξ1五1+…+ξ秒第页+…+ξ米五米=ξ1五1+…+ξ秒γ第页+2五第页+2+…+γ米五米+…+ξ米五米,

哪里λ我≥0且ξ我1≤时≥0我≤米现在将五秒给予λ秒=ξ秒γ秒，从那以后γ秒<0表示λ秒=ξ秒=0，所以问∩问(秒)也包含在维的子空间中米−1.

©2003英国皇家统计学会
本文根据牛津大学出版社标准期刊出版模式的条款出版和发行(https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
问题部分：
第条
下载所有幻灯片