2
$\begingroup$

Let $\phi(\cdot)$ and $\Phi(\cdot)$ be the probability and cumulative density functions, respectively, of a random variable with distribution $\text{N}(0,\,1)$. That is,
$$\forall x\in\mathbb{R}:\,\phi(x)=\frac{1}{\sqrt{2\cdot\pi}}\cdot\text{e}^{-x^2/2}$$ and
$$\forall x\in\mathbb{R}:\,\Phi(x)=\int_{-\infty}^{x}\phi(u)\,\text{d}u.$$ I was wondering if you could help me to compute
$$\int_{-\infty}^{w}\phi(x)\cdot\Phi(a+b\cdot x)\,\text{d}x\mbox{,}$$ where $(a,\,b,\,w)\in\mathbb{R}^3$ and $b\neq 0$, please.
Thanks a lot for your help.

$\endgroup$
3
  • 1
    $\begingroup$ How did this integral come up? Do you expect it to have some kind of closed-form solution? $\endgroup$
    – David Roberts
    Oct 20, 2017 at 6:42
  • $\begingroup$ Hi, @DavidRoberts. This integral is part of an algorithm I am programming in R-software. My purpose is to avoid the use of integrals in R-software because they reduce the efficiency of my algorithm. For this reason, I need to find a closed-form solution for this integral. $\endgroup$ Oct 20, 2017 at 10:09
  • $\begingroup$ Both Maple and Mathematica fail with it. Wanting is not harmful, harmful not want to. $\endgroup$
    – user64494
    Oct 20, 2017 at 10:51

2 Answers 2

3
$\begingroup$

This integral can be found in D. B. Owen (1980) A table of normal integrals, Communications in Statistics - Simulation and Computation, 9:4, 389-419: enter image description here

BvN denotes the bivariate normal probability function.

Check in R:

> a <- 2
> b <- 3
> w <- 5
> f <- function(x) dnorm(x)*pnorm(a+b*x)
> integrate(f, lower=-Inf, upper=w)
0.7364551 with absolute error < 1.3e-06
> 
> rho <- -b/sqrt(1+b^2)
> Sigma <- cbind(c(1,rho),c(rho,1))
> mvtnorm::pmvnorm(upper=c(a/sqrt(1+b^2), w), sigma=Sigma)
[1] 0.7364551
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"

Alternatively, you can express this integral with the Owen $T$-function:

> library(OwenQ)
> 1/2*(pnorm(a/sqrt(1+b^2))  + pnorm(w) - 2*OwenT(w, (b*w+a)/w) - 2*OwenT(-a/sqrt(1+b^2), (a*b+w*(1+b^2))/a) - (a <= 0))
[1] 0.7364551

Benchmark:

> library(mvtnorm)
> library(OwenQ)
> library(microbenchmark)
> 
> a <- 2
> b <- 3
> w <- 1
> 
> microbenchmark(
+   integral = integrate(function(x) dnorm(x)*pnorm(a+b*x), lower=-Inf, upper=w),
+   mvtnorm = {rho <- -b/sqrt(1+b^2); pmvnorm(upper=c(a/sqrt(1+b^2), w), sigma=cbind(c(1,rho),c(rho,1)))},
+   OwenT = 1/2*(pnorm(a/sqrt(1+b^2))  + pnorm(w) - 2*OwenT(w, (b*w+a)/w) - 2*OwenT(-a/sqrt(1+b^2), (a*b+w*(1+b^2))/a) - (a <= 0))
+ )
Unit: microseconds
     expr     min       lq      mean   median       uq      max neval cld
 integral  80.677  83.5860 116.97275  90.0240  93.0625 2878.062   100  b 
  mvtnorm 320.550 327.0625 339.22625 330.3975 336.0315  595.829   100   c
    OwenT  22.682  24.6360  28.89006  29.2685  31.9955   51.015   100 a  
$\endgroup$
4
  • 1
    $\begingroup$ @Student1981 I've just added a speed comparison in my answer. OwenT is the way to go. $\endgroup$ Oct 21, 2017 at 10:18
  • $\begingroup$ These are efficient reductions of the integral to standard libraries. But I'd say a solution with the bivariate normal CDF or OwenT is not in closed form; I would restrict the term "closed-form" to quantities that can be calculated with at most a single integral of elementary functions. $\endgroup$
    – user44143
    Oct 22, 2017 at 1:03
  • $\begingroup$ I appreciate your observation, @MattF. I believe I did not choose well my words when I wrote "closed-form solution". I will be more careful next time. However, for my purposes, the approach of Stéphane Laurent's answer is good enough. $\endgroup$ Oct 22, 2017 at 2:21
  • $\begingroup$ FYI, the condition in the OwenT derivation is wrong. It should be -(a/w<0) and not -(a<=0). The conclusion can be derived from 10,010.3 of Owen's paper and using the property $T(u,v) + T(uv,1/v) = \frac{1}{2} (\Phi(u) + \Phi(uv)) - \Phi(u) \Phi(uv) - \frac{1}{2} [v<0]$. $\endgroup$
    – jvdillon
    Feb 4, 2021 at 5:34
0
$\begingroup$

There is no known closed-form solution for this integral, though many have asked for it. However, for the case of $w=\infty$ you can use Geller and Ng (http://nvlpubs.nist.gov/nistpubs/jres/73B/jresv73Bn1p1_A1b.pdf), integral 4.3.13.

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.