Help me understand why
$\int_{-\infty}^{\infty}\frac{1}{2}[1+\operatorname{erf}(\frac{\theta-x}{\sqrt{2q^2}})]\frac{1}{\sqrt{2\pi\sigma^2}}{\exp(-\frac{(x-\mu)^2}{2\sigma^2})}dx \approx \frac{1}{2}[1+\operatorname{erf}(\frac{\theta-\mu}{\sqrt{2(q^2+\sigma^2)}})]$
This transformation used by Mark E. Glickma in "Parameter estimation in large dynamic paired comparison experiments" 1, but i can't find why it is so.
This is too long to be a comment.
Let $X$ and $Y$ be independent ${\rm N}(\mu,\sigma_2)$ and ${\rm N}(0,q^2)$ rv's, respectively.
Since $X+Y \sim {\rm N}(\mu,q^2+\sigma^2)$, it is equal in distribution to $Z + \mu$, where $Z \sim {\rm N}(0,q^2+\sigma^2)$. Hence,
$$
{\rm P}(X + Y \le \theta ) = {\rm P}(Z \le \theta - \mu ) = \frac{1}{{\sqrt {2\pi (\sigma ^2 + q^2 )} }}\int_{ - \infty }^{\theta - \mu } {e^{ - z^2 /[2(\sigma ^2 + q^2 )]} \,{\rm d}z} .
$$
On the other hand, by the law of total probability (conditioning on $X$), we have
$$
{\rm P}(X + Y \le \theta ) = \int_{ - \infty }^\infty {{\rm P}(Y \le \theta - x)\frac{1}{{\sqrt {2\pi \sigma ^2 } }}e^{ - (x - \mu )^2 /(2\sigma ^2 )} \,{\rm d}x}.
$$
Therefore,
$$
{\rm P}(X + Y \le \theta ) = \int_{ - \infty }^\infty {\bigg[\int_{ - \infty }^{\theta - x} {\frac{1}{{\sqrt {2\pi q^2 } }}e^{ - y^2 /(2q^2 )} \,{\rm d}y\bigg]\frac{1}{{\sqrt {2\pi \sigma ^2 } }}e^{ - (x - \mu )^2 /(2\sigma ^2 )} \,{\rm d}x} }.
$$
So,
$$
\int_{ - \infty }^\infty {\bigg[\int_{ - \infty }^{\theta - x} {\frac{1}{{\sqrt {2\pi q^2 } }}e^{ - y^2 /(2q^2 )} \,{\rm d}y\bigg]\frac{1}{{\sqrt {2\pi \sigma ^2 } }}e^{ - (x - \mu )^2 /(2\sigma ^2 )} \,{\rm d}x} },
$$
which may correspond to the left-hand side expression in the question, is equal to
$$
\frac{1}{{\sqrt {2\pi (\sigma ^2 + q^2 )} }}\int_{ - \infty }^{\theta - \mu } {e^{ - z^2 /[2(\sigma ^2 + q^2 )]} \,{\rm d}z},
$$
which may correspond to the right-hand side expression in the question (where $s^2$ should be $\sigma^2$).
This looks like a Laplace approximation. If $\sigma^2$ is small, the integral of $f(x)$ times a normal density is essentially $f(\mu)$. You could see this by expanding $f(x)$ in a power series centered at $\mu$ and keeping only the first two terms. The constant term gives the Laplace approximation and the second term integrates to zero by symmetry.
The Laplace approximation would give $$ \frac{1}{2} \left[1 + \textrm{erf}\left( \frac{\theta - \mu}{\sqrt{2q^2}} \right)\right] $$ which isn't quite the approximation in your question. But if $\sigma^2$ is sufficiently small, the two expressions are approximately equal. Perhaps the approximation in your question is a refinement of the Laplace approximation.
I believe this is in fact an exact result, as one can show by taking the derivative of the error function with respect to the additive constant and then reintegrating. See here for more details:
http://blitiri.blogspot.com/2012/11/gaussian-integral-of-error-function.html
with the substitutions: $y = x - \mu$ , $\alpha = 1/(2 \sigma ^2)$ , $\beta = (2 q^2)^{-1/2} \quad$ and $\gamma = \frac{\theta - \mu}{\sqrt{2 q^2}}$.
