Search code examples
rnormal-distributionnumerical-integration

Equivalent integrals not giving same numerical result with single and double integration

I have the following equivalent integrals involving the pdf and cdf of a standard normal distribution, denoted by $\phi(x)$ and $\Phi(x)$, respectively. For background, the integrals represent the probability distribution of the sample range of three observations, which are independent and identically distributed as a standard normal.

range_distribution

I code up the first double integral and the last single integral using the pracma package in R:

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

The output are:

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

Any ideas why these are giving different results when the underlying integrals are the same?

Solution

  • Your integrals are sensible to the choice of the bounds you take to replace infinity. In integral, you can use -Inf and Inf. In integral2, you can't. I would use the cubature package instead. When an integration bound is infinite, it automatically performs a change of variables to come down to a finite bound.

    library(cubature)

f <- function(xy) {
  x <- xy[1]; y <- xy[2]
  6 * dnorm(x) * dnorm(x+y) * (pnorm(x+y) - pnorm(x))
}

hcubature(
  f, lowerLimit = c(-Inf, 0), upperLimit = c(Inf, 0.2)
)

# 0.01096558