Closed form of the integral of the difference of two Gaussian CDFs?

Question

Problem

I'm trying to find the simplest form of the difference of two Guassian CDFs, i.e.

$$ \int_{-\infty}^\infty \left( \Phi\left(ax+b \right) - \Phi\left(cx+d \right) \right) dx $$

for $\Phi(\cdot)$ : standard norm CDF. Unfortunately there is no spoon-feeding solution for this problem in the Wiki page for the List of Integrals of Gaussian functions. For those who may flag this post, this post does not seem to address my question due to some misunderstandings on OP's part.

---

Try

But I tried with the indefinite integral of $\int \Phi(ax+b)dx$, which is given by

$$ \int \Phi(ax+b)dx = \frac{1}{b} \left ((a+bx)\Phi(a+bx) + \phi(a+bx)\right) + C $$

whose source is the above Wiki page. Thus

$$ \begin{aligned} \int \left( \Phi(ax+b) - \Phi(cx+d)\right) dx &= \frac{1}{b} \left ((a+bx)\Phi(a+bx) + \phi(a+bx)\right) \\ &- \frac{1}{d} \left ((c+dx)\Phi(c+dx) + \phi(c+dx)\right) + \tilde{C} \end{aligned} $$

but in view of definite integral, the tricky part is

$$ x \left[ \Phi(a+bx) - \Phi(c+dx) \right] \Large|\normalsize_{-\infty}^\infty $$

so my struggling does not seem to end up with something.

---

Speculation

In search of an answer, I kind of found out that if $a=c$,

$$ \int_{-\infty}^\infty \left( \Phi\left(ax+b \right) - \Phi\left(ax+d \right) \right) dx = (b-d)/a $$

without any proof. The way I found out this is via R. The following code illustrates how I speculated it(for some arbitrary mu, sd).

# mu, sd can be arbitrary numbers

a = 0.5; c = 0.5
b = 2; d = -4

lb = -1000; ub = 1000 # suff. large lower/upper bounds for test integration

f = function(theta){
  CDF1 = pnorm(a*theta + b, mean=mu, sd=sd)
  CDF2 = pnorm(c*theta + d, mean=mu, sd=sd)
  return(CDF1-CDF2)
}

f = Vectorize(f)
(val = integrate(f, lower=lb, upper=ub)$value)
# the result is (b+d)/a, for any b,d if a=c

But for the cases in which $a$, $c$ are different, I currently have no idea.

---

Question

I would like to find the general form for the cases in which $a \neq c$ with a proof. I may be missing a trivial part, but I would like to ask for a help. Any helping hands will be greeted. Thanks.

Answer 1

The $(b-d)/a$ result is correct when $a \gt 0.$ This post explains why. It generalizes the question broadly in order to reveal the underlying ideas.

---

Because $\Phi$ is not special in this regard, let's consider any distribution function $F_X$ for a random variable $X.$ Recall that by definition, $F_X(x) = \Pr(X\le x)$ for any real number $x.$

Suppose $X$ has a finite expectation $E_X.$ One expression for the expectation is

$$E[X] = \int_{0}^\infty (1 - F(x))\mathrm{d}x - \int_{-\infty}^0 F(x) \mathrm{d}x = \int_0^\infty (1 - (F(x) + F(-x)))\mathrm{d}x.\tag{*}$$

Let's study how this behaves under affine transformations of $X$:

1. $F_X(x+b) = \Pr(X \le x+b) = \Pr(X-b \le x) = F_{X-b}(x).$

2. When $a \gt 0,$ $F_X(ax) = \Pr(X \le ax) = \Pr(X/a \le x) = F_{X/a}(x).$

Thus, for positive $a$ and $c,$

$$\frac{E[X] - b}{a} = E[(X-b)/a] = \int_{0}^\infty (1 - F(ax+b))\mathrm{d}x - \int_{-\infty}^0 F(ax+b) \mathrm{d}x$$

and

$$\frac{E[X] - d}{c} = E[(X-d)/c] = \int_{0}^\infty (1 - F(cx+d))\mathrm{d}x - \int_{-\infty}^0 F(cx+d) \mathrm{d}x.$$

Subtracting the first from the second yields

$$\eqalign{\frac{b-E[X]}{a} - \frac{d - E[X]}{c} &=\int_{0}^\infty (1 - F(cx+d))\mathrm{d}x - \int_{-\infty}^0 F(cx+d) \mathrm{d}x \\&- \left(\int_{0}^\infty (1 - F(ax+b))\mathrm{d}x - \int_{-\infty}^0 F(ax+b) \mathrm{d}x\right) \\ &= \int_{0}^\infty (F(ax+b)-F(cx+d))\mathrm{d}x \\ &+ \int_{-\infty}^0 (F(ax+b) - F(cx+d)) \mathrm{d}x\\ &= \int_{-\infty}^\infty (F(ax+b)-F(cx+d))\mathrm{d}x .}$$

When $a \lt 0,$ replace $a$ by $-a = |a|$ and apply all results to the distribution of $-X.$

We have thereby established the following general result:

When $|E[X]| \lt \infty$ and $ac \ne 0,$ then $$\int_\mathbb{R} (F(ax+b)-F(cx+d))\,\mathrm{d}x = \frac{b-E[X]}{|a|} - \frac{d - E[X]}{|c|}.$$

In the question with positive $a$ and $c$ and $F=\Phi,$ we have $E[X] = 0,$ reducing the integral to $b/a - d/c.$ When $a=c$ this simplifies to $(b-d)/a,$ exactly as suggested in the question.

This result isn't quite the most general one: when $a=c,$ the result holds in the form $(b-d)/a$ even when $X$ does not have a finite expectation. This is most easily seen by integrating the quantile function $F^{-1}:$ see https://stats.stackexchange.com/a/18439/919.

---

In case any of these manipulations appear doubtful, here is numerical confirmation using a host of different distributions (some, like the Pareto and Student t, have infinite variance; others--the versions of a Binomial and Poisson distribution--are discrete). Each "Example ..." column corresponds to these randomly-chosen $(a,b,c,d):$

  Example 1 Example 2 Example 3
a    0.6267    0.8831   -0.3398
b   -0.7173   -0.4401   -0.5836
c   -0.9224    0.9378   -0.2596
d    1.0414   -0.3053   -0.7139

In Example 1 the signs of $a$ and $c$ differ; in Example 2 they are both positive; and in Example 3 they are both negative.

The output is

   Example 1 Example 2 Example 3   Method Distribution
1     -2.401   -0.1893    1.2598 Integral        Gamma
2     -2.401   -0.1893    1.2598  Formula        Gamma
3     -2.529   -0.2058    1.4871 Integral      Uniform
4     -2.529   -0.2058    1.4871  Formula      Uniform
5     -2.727   -0.2313    1.8382 Integral      Weibull
6     -2.727   -0.2313    1.8382  Formula      Weibull
7     -7.900   -0.8996   11.0329 Integral       Pareto
8     -7.900   -0.8996   11.0329  Formula       Pareto
9     -2.274   -0.1728    1.0326 Integral       Normal
10    -2.274   -0.1728    1.0326  Formula       Normal
11    -3.117   -0.2817    2.5314 Integral    Lognormal
12    -3.117   -0.2817    2.5314  Formula    Lognormal
13    -2.274   -0.1728    1.0326 Integral    Student t
14    -2.274   -0.1728    1.0326  Formula    Student t
15    -1.933   -0.1288    0.4265 Integral     Binomial
16    -1.933   -0.1287    0.4265  Formula     Binomial
17    -2.444   -0.1948    1.3356 Integral      Poisson
18    -2.444   -0.1948    1.3356  Formula      Poisson

Each pair of lines shows the integral's value followed by the formula's value; they agree in every case.

Here is the code that performed these computations. Notice how the expectation $E[X]$ is carried out with the integral $(*)$ in the function g. The blind integration of the discrete distribution functions in g can be a little delicate; this is handled by increasing the default number of subdivisions from 100 to 1000, but could be further improved by using a finite lower limit of integration (thereby giving the routine a decent hint concerning the scale of the calculation).

#
# Compute the original integral numerically.
#
g <- function(a,b,c,d, F.=pnorm, ...) {
  integrate(function(x) F.(a*x + b) - F.(c*x + d), -Inf, Inf, ...)$value
}
#
# Apply the formula.  This requires knowing or finding E_F, the expectation
# of `F.`.  Here we find that expectation with a numerical integration.
#
g. <- function(a,b,c,d, F.=pnorm) {
  expectation <- integrate(function(x) 1 -  (F.(-x) + F.(x)), -Inf, 0)$value
  (1/abs(c) - 1/abs(a)) * expectation + b/abs(a) - d/abs(c)
}
#
# This is the Pareto CDF.
#
pPareto <- function(x, alpha, x.min) ifelse(x <= x.min, 0, 1 - (x/x.min)^(-alpha))
#
# Test a bunch of different distributions.
#
distributions <- list(Gamma = function(x) pgamma(x, 0.25),
                      Uniform = punif,
                      Weibull = function(x) pweibull(x, 2, 1),
                      Pareto = function(x) pPareto(x, 1.1, 1),
                      Normal = pnorm,
                      Lognormal = plnorm,
                      `Student t` = function(x) pt(x, 1.1),
                      Binomial = function(x) pbinom(x+4, 10, 1/3),
                      Poisson = function(x) ppois(x, 1/3))
#
# Create some random sets of (a,b,c,d) values.
#
set.seed(17)
coeffnames <- c("a","b","c","d")
args <- lapply(1:3, function(i) {x <- as.list(rexp(4)-1); names(x) <- coeffnames; x})
names(args) <- paste("Example", seq_along(args))
print(matrix(unlist(args), 4, dimnames=list(coeffnames, names(args))), digits=4)

#
# Conduct the tests.
#
Results <- do.call(rbind, lapply(names(distributions), function(s) {
  G <- distributions[[s]]
  X <- as.data.frame(rbind(sapply(args, do.call, 
                                  what=function(...) g(..., F.=G, subdivisions=1000L)),
       sapply(args, do.call, what=function(...) g.(..., F.=G))))
  X$Method <- c("Integral", "Formula")
  X$Distribution <- s
  X
}))
print(Results, digits=4)

Answer 2

A geometrical intuition

A geometrical intuition to accompany Whuber's answer is the following:

It relates to answers to other questions here and here.

The mean of the variable relates to the area of the gray striped surfaces, which can be computed in two directions

• The vertical stripes: as an integral of the quantile function with the quantiles $dp$ as integrand

  $$E[X] = \int_0^1 Q_X(p) dp$$

• The horizontal stripes: as an integral with the CDF and the variable $dz$ as integrand

  $$E[X] = -\int_{-\infty}^0 F_X(z) dz + \int_{0}^{\infty} 1-F_X(z) dz$$

Then the difference $E[X]-E[Y]$ is

$$\begin{array}{} E[X]-E[Y] &=& \int_{0}^1 Q_X(p)-Q_Y(p) dp \\ \end{array}$$

which relates to the difference of the quantile functions, and the blue and red areas marked by the vertical stripes.

But we can just as well compute this in the horizontal direction by taking the difference between the CDF functions.

$$\begin{array}{} E[X]-E[Y] &=& -\int_{-\infty}^0 F_X(z)-F_Y(z) dz + \int_{0}^{\infty} -F_X(z)+F_Y(z) dz \\ &=& - \int_{-\infty}^\infty F_X(z)-F_Y(z) dz \end{array}$$

Applied to Gaussian CDF

With the general equation

$$\begin{array}{} E[X]-E[Y] &=& - \int_{-\infty}^\infty F_X(z)-F_Y(z) dz \end{array}$$

you will get for the CDF's of the Gaussian distribution $\Phi\left(ax+b \right)$ which has the mean $\mu(a,b) = -b/a$

$$\int_{-\infty}^\infty \left( \Phi\left(ax+b \right) - \Phi\left(cx+d \right) \right) dx = -[\mu(a,b) - \mu(c,d)] = -b/a - c/d$$