How can I compute $\int F(x \mid a,b)f(x \mid w,z) {}dx$ in closed form? Suppose $F$ is the cumulative distribution function of the normal distribution with mean $a$ and standard deviation $b$, and suppose $f$ is the probability density function of the normal distribution with mean $w$ and standard deviation $z$. How can I calculate this integral in closed form?
$$\int F(x \mid a,b)f(x \mid w,z) {}dx$$
It seems like my question is very similar to this one, but crucially I'm not using that standard versions of the cdf or pdf. I tried converting my functions to their standard forms, but ended up with this:
$$\int \Phi\left(\frac{x-a}{b}\right)\frac{\phi\left(\frac{x-w}{z}\right)}{z} {}dx$$
And now I'm really confused!
P.S.: If the above does have a closed-form simplification, my true problem is actually closer to this:
$$\int F(2a - \xi - x \mid a,b)f(x \mid w,z) {}dx$$
 A: $$\int F(x \mid a,b)f(x \mid w,z) {}dx=E[F(Y \mid a,b)]=\Pr(X<Y)$$ 
where $Y \sim F(\cdot \mid w,z)$ is independent of $X \sim F(\cdot \mid a,b).$
Now, $X-Y \sim {\cal N}(a-w, z^2+b^2)$. The result is $\boxed{F(0 \mid a-w, \sqrt{z^2+b^2})}$.
Just to check:
> a <- 1; b <- 1
> w <- 2; z <- 2
> integrate(function(x) pnorm(x,a,b)*dnorm(x,w,z), lower=-Inf, upper=Inf)
0.6726396 with absolute error < 2.2e-09
> pnorm(0, a-w, sqrt(b^2+z^2))
[1] 0.6726396

Your last question directly follows from the first one by noting that
$$
F(2a - \xi - x \mid a,b) = 1 - F(x \mid 2a-\xi-a,b),
$$
therefore the result for your last integral is $\boxed{1-F(0 \mid a-\xi-w, \sqrt{z^2+b^2})}$.
> ksi <- 3
> integrate(function(x) pnorm(2*a-ksi-x,a,b)*dnorm(x,w,z), lower=-Inf, upper=Inf)
0.03681913 with absolute error < 0.00012
> 1- pnorm(0, a-ksi-w, sqrt(b^2+z^2))
[1] 0.03681914

A: I'll gibe you the special case, which may point you to the general case, if it exists. 
If the parameters of your distributions were the same, then we have
$$\int F(x|\mu,\sigma)f(x|\mu,\sigma)dx=\int F(x|\mu,\sigma)dF(x|\mu,\sigma)=F^2(x|\mu,\sigma)/2+Const$$
When parameters are not the same, then it gets hairy. There's got to be a way to convert them to standard normal, then it's be the same parameters (0,1). I think it'll give you solution with constants and scaling
So, for different parameters, something like this should be true:
$$F(x|\mu_2,\sigma_2)=   F(\frac{x-\mu_2}{\sigma_2}\sigma+\mu|\mu,\sigma)$$
Then the integral becomes:
$$\int F(x|\mu_2,\sigma_2)f(x|\mu,\sigma)dx=\int F(\frac{x-\mu_2}{\sigma_2}\sigma+\mu|\mu,\sigma)dF(x|\mu,\sigma)$$
or
$$=\int y(a x +b) dy(x)=?$$
This got to have some beautiful solution, Stiltjes integral or something
UPDATE: I vaguely remember that this type of integrals was used a lot on game theory. The stochastic dominance of the second kind and stuff like that is all based on these integrals
