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