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$ – theHigherGeometer Oct 20 '17 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$ – Student1981 Oct 20 '17 at 10:09
  • $\begingroup$ Both Maple and Mathematica fail with it. Wanting is not harmful, harmful not want to. $\endgroup$ – user64494 Oct 20 '17 at 10:51
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$ – Stéphane Laurent Oct 21 '17 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$ – Matt F. Oct 22 '17 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$ – Student1981 Oct 22 '17 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 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, privacy policy and cookie policy

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