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!
 A: 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.
A: $$\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).
A: 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. 
A: 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
  #      - Phi(x) Phi(a/rho)                    # wikipedia page on Owen's T.
  #      - [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
  # logic to preserve gradients.
  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))
x = tf.boolean_mask(x,(x>lo)&(x<hi))
[
  tf.reduce_mean(tfd.Normal(0,1).cdf(x)) * (d.cdf(hi) - d.cdf(lo)),
  gaussian_expected_standard_gaussian_cdf(m,s,lo,hi)
]

