$$X \sim \mathcal N(\mu, \sigma^2)$$ $$Y \sim Bern(p)$$
$$Z = XY$$
I then have that the pdf of $Z$ is:
$$f(z) = (p-1)\delta(z) + p g(z)$$
where $g(z)$ is the pdf of the normal distribution.
I then need to find the sum of $Q$ of these distributions and finally have to find the expected value of the absolute value of the final resulting distribution.
I do not know where to go from here. I have tried doing it with MGFs but I hit a brick wall with divergent integrals. Any help please?
EDIT:
I forgot to add that $\frac{Q}{2}$ of these distributions will have the normal with $-\mu$ and the other half with $\mu$. I have tried doing this with a Monte Carlo simulation and then with a similar method to what Henry suggested in the comments.
My simulation has the following results:
| Q | mean | stdev |
|---|---|---|
| 2 | 3.406746 | 3.146736 |
| 4 | 5.156385 | 4.041445 |
| 6 | 6.393295 | 4.865093 |
| 8 | 7.391076 | 5.586568 |
| 10 | 8.283423 | 6.254462 |
| 12 | 9.075857 | 6.846168 |
| 14 | 9.788690 | 7.383987 |
| 16 | 10.471054 | 7.905363 |
| 18 | 11.102388 | 8.384783 |
| 20 | 11.715745 | 8.837827 |
And the equation I got to with the various combinations of folded normal distributions is, where I have used that $p = 0.5$ (if anyone has a more generic solution even better):
FURTHER EDIT
I got it for $p = 0.5$, I had made a silly mistake at the start. Thank you all.
Denote $f(a)=E(|Z-a|)$ if $Z\sim N(0,1)$ Denoting $\Phi(a)=\Pr(Z<a)$ we get $$f(a)=2a\Phi (a)-a+\sqrt{\frac{2}{\pi}}e^{-a^2/2}.$$ Now let $X_1,\ldots,X_n$ and $Y_1,\ldots, Y_n$ independent such that $X_i$ is Bernoulli of mean $p$ and $Y_i=\mu+\sigma Z_i$ with $Z_i\sim N(0,1). $ Note that $$|Y_1+\cdots+Y_k|=|k\mu+\sigma (Z_1+\cdots+Z_k)|\sim \sqrt{k}\sigma\left|Z_1+\frac{\sqrt{k}}{\sigma }\mu\right|$$ since $Z_1+\cdots+Z_k\sim N(0,k)$ and therefore $E(|Y_1+\cdots+Y_k|)=\sqrt{k}\sigma f(-\frac{\sqrt{k}}{\sigma }\mu)=g_k.$ To conclude$$E(|X_1Y_1+\cdots+X_nY_n|)=\sum_{k=0}^nC^k_np^k(1-p)^{n-k}g_k.$$