# Integrating $\frac1{2\sigma\sqrt{2\pi}}\int_{-\infty}^\infty(1+erf(\frac{z-\mu_i}{\sigma_i\sqrt2}))\exp( -\frac12(\frac{z-\mu}\sigma)^2)dz$

I am interested in the integral of $$\frac{1}{2\sigma \sqrt{2\pi}}\int_{-\infty}^\infty \left( 1+\operatorname{erf}\left(\frac{z-\mu_i}{\sigma_i\sqrt{2}}\right)\right)\exp\left( -\frac{1}{2}\left( \frac{z-\mu}{\sigma}\right)^2\right)dz$$

where $$\operatorname{erf}$$ is the error function defined as the CDF of the standard normal (up to a constant).

In this case $$\sigma, \mu, \mu_i, \sigma_i$$ are constants.

I truly have no clue how to even start something like this. Is Mathematica even capable of this type of integral?

For those curious, this is evaluated when considering $$Pr(f_i \leq f \leq f_j)$$ where $$f_j, f_i, f$$ follow a normal distribution with different parameters. The full integral is

$$\frac{1}{2\sigma \sqrt{2\pi}}\int_{-\infty}^\infty \left( 1+erf\left(\frac{z-\mu_i}{\sigma_i\sqrt{2}}\right)\right)\exp\left( -\frac{1}{2}\left( \frac{z-\mu}{\sigma}\right)^2\right)dz \\ - \frac{1}{4 \sigma \sqrt{2 \pi}} \int_{-\infty}^\infty \left( 1+erf \left(\frac{z-\mu_i}{\sigma_i \sqrt{2}}\right)\right)\left(1+erf\left(\frac{z-\mu_j}{\sigma_j \sqrt{2}}\right)\right)\exp\left(-\frac{1}{2}\left(\frac{z-\mu}{\sigma}\right)^2\right)dz$$

But I thought that might be inappropriate ;)

• You are looking for the closed-form of $\int_{\mathbb R} e^{-(x-x_0)^2}\operatorname{erf}x\,dx$, which might not exist generally. Feb 17 '21 at 6:09
• @Julian This is a gentle reminder to consider accepting an answer if your problem has been resolved.
– Sal
Feb 23 '21 at 20:25

Integral of one error function and Gaussian

Mathematica does not evaluate the integral. Surprisingly, there is an exact result. By changing variables, $$z-\mu=x$$, your integral of $$\operatorname{erf} \times \text{Gaussian}$$ may be written as

$$I(b)=\int\limits_{-\infty}^\infty dx \ \operatorname{erf}(ax+b) \exp \left(-cx^2 \right)$$

Since $$\partial_t \operatorname{erf}(t)= \frac{2}{\sqrt{\pi}}e^{-t^2}$$, differentiating with respect to $$b$$ produces

$$\partial_b I(b)=\int\limits_{-\infty}^\infty dx \ \frac{2}{\sqrt{\pi}} \exp\left( -(ax+b)^2-cx^2 \right)$$

This is just a Gaussian integral over $$x$$, the result is

$$\partial_b I(b)= \frac{2}{\sqrt{a^2+c}} \exp \left( -\frac{b^2 c}{a^2 +c} \right)$$

Now we integrate to recover $$I(b)$$

$$I(b)-I(0)=\int\limits_0^b \ dt \frac{2}{\sqrt{a^2+c}} \exp \left( - t^2\frac{ c}{a^2 +c} \right)$$

$$I(0)=\int dx \ ( \text{odd function})=0$$, and the integral on the right is simply the definition of $$\operatorname{erf}$$, so we have

$$I(b)=\sqrt{\frac{\pi}{c}} \operatorname{erf} \left( b \sqrt{\frac{c}{a^2+c}} \right)$$

$$\tag*{\blacksquare}$$

Integral of two error functions and Gaussian

We may use a similar method to find that

$$\int\limits_{-\infty}^\infty dx \ \operatorname{erf}(a_1 x) \operatorname{erf}(a_2 x) e^{-cx^2} = \frac{2}{\sqrt{c\pi}} \tan^{-1}\left( \frac{a_1 a_2}{ \sqrt{c(a_1^2+a_2^2+c)}} \right)$$

However, once we introduce $$b\neq 0$$, we meet the more complicated integral in an intermediate step

$$\int\limits_0^b dt \ \operatorname{erf}(\alpha t + \beta) e^{-\gamma t^2}$$

Update! Using the result from this post, we can make progress. Let

$$J(\beta)=\int\limits_{-\infty}^\infty dx \ \operatorname{erfc}(\alpha x+\beta) \operatorname{erf}(a x+b) \exp \left(-cx^2 \right)$$

Then $$J(\infty)=0$$. Differentiating

$$\partial_{\beta} J(\beta)=-\frac{2}{\sqrt{\pi}} \int\limits_{-\infty}^\infty dx \ \operatorname{erf}(a x+b) \exp \left(-cx^2 - (\alpha x +\beta)^2 \right)$$

Which we can evaluate using our first result

$$\partial_{\beta} J(\beta)=-\frac{2}{\sqrt{\pi}} \sqrt{\frac{\pi}{c+\alpha^2}} \operatorname{erf}\left( \left(b-\frac{a \alpha \beta}{c+\alpha^2} \right) \sqrt{\frac{c+\alpha^2}{c+\alpha^2+a^2}} \right) \exp \left( -\beta^2 \left( 1+\frac{\alpha^2}{c+\alpha^2}\right) \right)$$

With $$\Lambda=\frac{2}{\sqrt{\pi}} \sqrt{\frac{\pi}{c+\alpha^2}}$$, $$A=-\frac{a \alpha }{c+\alpha^2} \sqrt{\frac{c+\alpha^2}{c+\alpha^2+a^2}}$$, $$B=b\sqrt{\frac{c+\alpha^2}{c+\alpha^2+a^2}}$$, and $$C=\left( 1-\frac{\alpha^2}{c+\alpha^2}\right)$$. Now we integrate

$$J(\infty)-J(\beta)=-\int\limits_\beta^\infty \ dt \Lambda \operatorname{erf}(At+B)e^{-Ct^2}$$

$$J(\beta)=\int\limits_\beta^\infty dt \ \Lambda \operatorname{erf}(At+B)e^{-Ct^2}= \frac{\Lambda}{\sqrt{2C}} \int\limits_{\sqrt{2C}\beta}^\infty du \ \operatorname{erf}(Au/\sqrt{2C}+B)e^{-u^2/2}$$

Which may be written in terms of the 'Generalized Owen-T' function

$$J(\beta)=\frac{2 \sqrt{\pi} \Lambda}{\sqrt{C}}T\left( \sqrt{2C}\beta, A/\sqrt{C},\sqrt{2}B\right)$$

Using identity ii we have

$$J(\beta)=\frac{2 \sqrt{\pi} \Lambda}{\sqrt{C}} \left[ \frac{1}{2\pi} \left( \tan^{-1}(A/\sqrt{C}) -\tan^{-1}(A/\sqrt{C} + B/\beta\sqrt{C}) \\ -\tan^{-1}(B^{-1}(\beta \sqrt{C} +AB/\sqrt{C} + A^2 \beta /\sqrt{C} )) \right) +\frac{1}{4} \operatorname{erf}(B (1+A^2/C)^{-1/2}) \\ +T(\beta \sqrt{2C},(A\beta^{-1}+B)/\beta \sqrt{C})+T(\sqrt{2}B (1+A^2/C)^{-1/2}, B^{-1}(\beta \sqrt{C} +AB/\sqrt{C} + A^2 \beta /\sqrt{C} ) \right]$$

Where $$T$$ is the Owen T function. I've checked it numerically, I hope there are no typos here.

Approximate integral of two or more error functions with Gaussian

We can approximate the integral using Laplace's method. I'll demonstrate with two error functions, but in principle you could have more. Let

$$I(c)=\int\limits_{-\infty}^\infty dx \ \operatorname{erf}(a_1 x+ b_1) \operatorname{erf}(a_2 x+ b_2) e^{-cx^2}$$

Expand the product of $$\operatorname{erf}$$s around the Gaussian peak: $$x=0$$

$$\operatorname{erf}(a_1 x+ b_1) \operatorname{erf}(a_2 x+ b_2) = \operatorname{erf}( b_1) \operatorname{erf}( b_2) + x \left( \text{junk} \right) + x^2 \frac{2 e^{-b_1^2-b_2^2}}{\pi} \left( 2 a_1 a_2 -a_2^2 b_2 e^{b_1^2} \sqrt{\pi} \operatorname{erf}(b_1) - a_1^2b_1 e^{b_2^2} \sqrt{\pi} \operatorname{erf}(b_2) \right) + \mathcal{O}(x^3)$$

Term by term integration produces the asymptotic series. The $$(\text{junk})$$ and all $$x^{\text{odd}}$$ terms integrate to zero. Noting that

$$\int\limits_{-\infty}^\infty dx \ x^{2n} e^{-cx^2} = c^{-n-1/2} \Gamma(n+1/2)$$

For fixed $$a_i, \ b_i$$ we have

$$I(c) \sim \sqrt{\frac{\pi}{c}} \operatorname{erf}( b_1) \operatorname{erf}( b_2) \\ +\frac{1}{c^{3/2}} \frac{ e^{-b_1^2-b_2^2}}{\sqrt{\pi}} \left( 2 a_1 a_2 -a_2^2 b_2 e^{b_1^2} \sqrt{\pi} \operatorname{erf}(b_1) - a_1^2b_1 e^{b_2^2} \sqrt{\pi} \operatorname{erf}(b_2) \right) \ \ , \ \ c \to \infty$$

Which is the two term approximation.

• Holy (explicative)! I am amazed :). Would a similar technique apply for the second term, where we have as the integrand $\text{erf}(ax+b)\text{erf}(cx+d)\exp(-gx^2)$ where $a,b,c,d,g$ are constants? Feb 19 '21 at 1:18
• Great answer @Sal! Many thanks! Jul 27 '21 at 11:17