5
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ I think it's also true that for any $a,c$ if $b=d$ then $\int_{-\infty}^{\infty}\Phi(ax+b)-\Phi(cx+d)\,dx=0$? $\endgroup$ – Ed P Jan 25 at 8:44
  • 2
    $\begingroup$ @Stats The integral can diverge, too, when $ac \le 0.$ $\endgroup$ – whuber Jan 25 at 14:24
  • 1
    $\begingroup$ FWIW, there is some confusion in this post concerning the roles of $a,b,c,d.$ I don't think that makes the question incomprehensible, but readers should be careful. $\endgroup$ – whuber Jan 25 at 16:31
5
$\begingroup$

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)
| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Thank you for a beautiful answer as always. $\endgroup$ – moreblue Jan 26 at 1:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.