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.
$\begingroup$
$\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 ♦Commented 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$– Student1981Commented Oct 20, 2017 at 10:09
-
$\begingroup$ Both Maple and Mathematica fail with it. Wanting is not harmful, harmful not want to. $\endgroup$– user64494Commented Oct 20, 2017 at 10:51
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
4
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:
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
answered Oct 20, 2017 at 16:12
-
1$\begingroup$ @Student1981 I've just added a speed comparison in my answer.
OwenTis the way to go. $\endgroup$ Commented 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$– user44143Commented 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$ Commented 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$– jvdillonCommented Feb 4, 2021 at 5:34
$\begingroup$
$\endgroup$
0
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.