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 ;)
 A: 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.
