等价正态分布积分单、双数值计算结果不一致原因排查

问题描述 #

我们有理论等价的积分，用于表示三个独立同分布标准正态样本的极差概率分布，涉及标准正态分布的概率密度函数φ(x)（对应R语言的dnorm()）和累积分布函数Φ(x)（对应R语言的pnorm()）。使用pracma包实现单积分和二重积分后，得到的结果却完全不同：

用户的实现代码：

library(pracma)

t = 0.2
single_int <- function(x) 3*dnorm(x)*(pnorm(x+t)-pnorm(x))^2
I1 = integral(single_int,-999,999)       

double_int <- function(x,y) 6*dnorm(x)*dnorm(x+y)*(pnorm(x+y)-pnorm(x))
I2 = integral2(double_int,-999,999,-999,t)$Q

运行结果：

> I1
[1] 0.01096555
> I2
[1] 0 

问题原因 #

1. 二重积分的区间设置错误 #

从理论等价性来看，二重积分中y的取值范围应该是0到t，而不是你设置的-999到t。当y < 0时，pnorm(x+y)-pnorm(x)为负数，这部分积分会和y > 0的部分相互抵消，最终导致整体积分结果为0。而单积分里的平方项消除了符号影响，所以能得到正确结果。

2. 过宽的积分区间导致数值问题 #

把积分下限设为-999虽然试图逼近负无穷，但数值积分中过宽的区间会增加计算负担甚至引入误差。标准正态分布的几乎全部概率都集中在[-4,4]区间内，用这个区间足以保证计算精度，同时提升效率。

修正后的代码 #

library(pracma)

t = 0.2
# 优化单积分区间
single_int <- function(x) 3*dnorm(x)*(pnorm(x+t)-pnorm(x))^2
I1 = integral(single_int, -4, 4)       

# 修正二重积分的y区间为0到t
double_int <- function(x,y) 6*dnorm(x)*dnorm(x+y)*(pnorm(x+y)-pnorm(x))
I2 = integral2(double_int, -4, 4, 0, t)$Q

修正后运行结果：

> I1
[1] 0.01096555
> I2
[1] 0.01096554

内容的提问来源于stack exchange，提问作者stats134711