Definite integral of a product of normal pdf and cdf

Denote the pdf of the standard normal distribution as $\phi(x)$ and cdf as $\Phi(x)$. Does anyone know how to calculate $\int_{-\infty}^y \phi(x)\Phi(\frac{x−b}{a})dx$?

Notice that this question is similar to an existing one,

https://mathoverflow.net/questions/101469/integration-of-the-product-of-pdf-cdf-of-normal-distribution

the only difference being that I'm computing the integral over $(-\infty, y)$ for some real $y$, rather than over the entire real line.

Thank you!

• 1) not a research level question 2) what does "ba" mean??? – Alexey Jun 30 '13 at 18:52
• Where does the problem come from? Why do you need this result? – Davide Giraudo Jun 30 '13 at 19:20
• Alexey: 1) It's difficult for me to judge the level and complexity of this question. It is likely that you are in a better position to do so, so I'll take your word for it. 2) This was a typo, thank you very much for pointing it out. I have corrected the question. The argument of $\Phi$ should read as $(\frac{x - b}{a})$, not $(x - ba)$ as it was. Here $a, b$ are real constants, $a \neq 0$. – user36247 Jul 1 '13 at 10:40
• Davide: This integral has appeared in several contexts in my investigation of Brownian motion over compact intervals, hence the need to integrate over $(-\infty, y)$. (In fact, I need to find this integral over $(y, z)$ for some real $y$ and $z$.) I don't know whether it has an analytic solution. Hopefully it does. I would very much appreciate your help. – user36247 Jul 1 '13 at 10:43
• I had a very similar question, in fact the same here: math.stackexchange.com/questions/2746095/… Sadly, the answer is that the definite integral does not have closed form solution. I offered also graph interpretation in L_2 space on the picture, to demonstrate another view on the problem. – Jan Muchna May 4 '18 at 7:06

As already explained, when $a\gt0$ the full integral is $1-\Phi\left(b/\sqrt{a^2+1}\right)$. The same approach shows that the integral considered here is $$I(y)=P(Y\leqslant(X-b)/a,X\leqslant y),$$ where $(X,Y)$ are i.i.d. standard normal, that is, $$I(y)=P(aY+b\leqslant X\leqslant y).$$ I see no reason to expect more explicit formulas.

• Hi, Did, I doubt there is a more explicit formula. Someone told me this integral may appear in an old paper by Owen. I'll look it up and let you know. – user36247 Jul 3 '13 at 16:06
• By all means, please do. :-) – Did Jul 3 '13 at 19:03

$$\int_{-\infty}^y \phi(x) \Phi(\frac{x-b}{a})dx = BvN\left[\frac{-b}{\sqrt{a^2+1}}, y; \rho= \frac{-1}{\sqrt{a^2+1}}\right]$$

where $BvN(w, z; \rho)$ is the bivariate normal cumulative with upper bounds $w$ and $z$, and correlation $\rho$.

For reference, see equation (10,010.1) in Owen (Comm. in Stat., 1980).

If you check out the integral tables in section 4.2 and 4.3 of http://nvlpubs.nist.gov/nistpubs/jres/73B/jresv73Bn1p1_A1b.pdf, you will find what you need to get this done. I used equations 4.3.2 and 4.2.1 (latter is same as eqn 7.4.36 in Abranovitz and Stegun). Just change variables on the error function and complete the square on the exponential. You will end up with one 4.3.2 integral, one 4.2.1, and one simple erf(x) integral of the exponential square.

I am a bioinformatician. This problem occured for me in the context of statistics. I was trying to compute conditional probabilities to input in my factor graph model.

I verified equation 4.2.1 from the source litterature. I will have to doublecheck if 4.3.2 is correct (via numerical integration), since this solution is original to this work.

• Hi Hughes. Could you document your solution here? (Ideally with derivation steps.) Thanks! – jvdillon Feb 2 at 16:36

What follows isn't exactly the answer to your question but should be what you need after you make the necessary change of variables.

\begin{align} &\int_{l}^{h} \textrm{d} x\, \frac{1}{\sigma}\phi\left(\frac{x-\mu}{\sigma}\right) \Phi(x)\nonumber\\ =& \tfrac{1}{2}(\Phi(z_h) - \Phi(z_l)) \\ & - \tfrac{1}{2}\left[\frac{\mu}{z_h} < 0\right] \\ & + \tfrac{1}{2}\left[\frac{\mu}{z_l} < 0\right] \\ & - T\left(z_h, \frac{h}{z_h}\right) \\ & + T\left(z_l, \frac{l}{z_l}\right) \\ & - T\left(\frac{\mu}{\rho}, \frac{\mu \sigma + z_h \rho^2}{\mu}\right) \\ & + T\left(\frac{\mu}{\rho}, \frac{\mu \sigma + z_l \rho^2}{\mu}\right) \end{align} where \begin{align} \rho &=\sqrt{1+\sigma^2} \\ z_l &= \frac{l - \mu}{\sigma} \\ z_h &= \frac{h - \mu}{\sigma} \\ \end{align} and where $$\phi,\Phi$$ are the pdf and cdf of the standard Normal distribution (i.e. Normal with mean $$0$$ and std dev $$1$$) and $$T$$ denotes the Owen's T function.

Here's an implementation of this integral using TensorFlow Probability; it should be easily adaptable to any numerical backend which offers Owen's T function.

import tensorflow as tf
import tensorflow_probability as tfp
from tensorflow_probability.python.internal import dtype_util
tfd = tfp.distributions

def gaussian_expected_standard_gaussian_cdf(m, s, lo, hi):
"""Computes:
E[ Normal(0,1).cdf(X) | lo<X<hi, X~Normal(m,s)] * (
Normal(m,s).cdf(hi) - Normal(m,s).cdf(lo))
"""
m, s, lo, hi, dtype = _prepare_args(m, s, lo, hi)
d = tfd.Normal(dtype(0), 1)
# We derive our implementation by beginning with the result 10,010.3 (pg 403) of "A Table of Normal Integrals" by DB Owen.
#
# rho = sqrt(1+b^2)
#
# int phi(x) Phi(a + bx) dx                   # Result #10,010.3, pg 403.
# =   T(x, a / (x rho)) + T(a / rho, x rho / a)
#   - T(x, (a+bx)/x)
#   - T(a / rho, (a b + x rho**2)/a)
#   + Phi(x) Phi(a / rho)
#   + C
# =   ([Phi(x) + Phi(a/rho)]/2                # Property 2.7, pg 414. See also
#      - [a/x < 0]/2)
#   - T(x, (a+bx)/x)
#   - T(a / rho, (a b + x rho**2)/a)
#   + Phi(x) Phi(a / rho)
#   + C
# =  Phi(x)/2 - [a/x < 0]/2
#   - T(x, (a+bx)/x)
#   - T(a / rho, (a b + x rho**2) / a)
#   + C + Phi(a/rho)/2
#
# Evaluate the indefinite integral from u to v:
#
#  + (Phi(v) - Phi(u))/2
#  - [a/v < 0]/2
#  + [a/u < 0]/2
#  - T(v, (a+bv)/v)
#  + T(u, (a+bu)/u)
#  - T(a / rho, (a b + v rho**2) / a)
#  + T(a / rho, (a b + u rho**2) / a)
#
# Make substitutions a=m, b=s and change of variables
# v=(hi-m)/s=:zhi  and u=(lo-m)/s:=zlo:
#
#  + (Phi(zhi) - Phi(zlo)) / 2
#  - [m / zhi < 0]/2
#  + [m / zlo < 0]/2
#  - T(zhi, hi / zhi)
#  + T(zlo, lo / zlo)
#  - T(m / rho, (m s + zhi rho**2) / m)
#  + T(m / rho, (m s + zlo rho**2) / m)
rho2 = s**2 + 1.
m_over_rho = m * tf.math.rsqrt(rho2)
zhi = (hi - m) / s
zlo = (lo - m) / s
return (
0.5 * (d.cdf(zhi) - d.cdf(zlo))
- tf.where(tf.equal(m, 0) | tf.equal(m < 0, zhi < 0),
dtype(0.), dtype(0.5))
+ tf.where(tf.equal(m, 0) | tf.equal(m < 0, zlo < 0),
dtype(0.), dtype(0.5))
- _safe_owens_t(zhi, hi, zhi)
+ _safe_owens_t(zlo, lo, zlo)
- _safe_owens_t(m_over_rho, m * s + zhi * rho2, m)
+ _safe_owens_t(m_over_rho, m * s + zlo * rho2, m)
)

# We need to be able to call Owen's T with some potentially "odd"
# inputs; the following function ensures no NaNs and preserves gradients
# using the "double where trick."
def _safe_owens_t(h, numer, denom):
# Property 2.4 (pg 414) states that:
#   T(h, inf) = (1 - Phi(|h|))/2
# We use this property to skip the Owen's T thus avoiding nans.
# We'll interpret 0/x=0 and inf/inf=inf; use what may seem like convoluted
is_numer_zero = tf.equal(numer, 0)
is_denom_nonzero = tf.not_equal(denom, 0)
has_finite = tf.math.is_finite(denom) | tf.math.is_finite(numer)
safe_denom = tf.where(is_denom_nonzero & has_finite, denom, 1.)
dtype = dtype_util.as_numpy_dtype(h.dtype)
return tf.where(
is_numer_zero | is_denom_nonzero,
tfp.math.owens_t(h, numer / safe_denom),
0.5 * tfd.Normal(dtype(0), 1).survival_function(tf.math.abs(h)))

def _prepare_args(m, s, lo, hi=None):
dtype = dtype_util.common_dtype([m, s], dtype_hint=tf.float32)
m = tf.convert_to_tensor(m, dtype, name='m')
s = tf.convert_to_tensor(s, dtype, name='s')
lo = tf.cast(lo, dtype, name='lo')
dtype = dtype_util.as_numpy_dtype(dtype)
if hi is None:
return m, s, lo, dtype
hi = tf.cast(hi, dtype, name='hi')
return m, s, lo, hi, dtype


You can test the code like:

m,s = -0.25,1.43
lo,hi = -np.inf,1.
d = tfd.Normal(m,s)
x = d.sample(int(100e6))