Symbolic integration of a indexed product term?

Question

I am trying to perform the symbolic integration of the following expression $$\int \prod_{i=1}^n p(x_i | z) p(z) dz$$ where $$p(x_i|z) = \frac{\exp\left(\tfrac{-(x-z)^2}{2\sigma_i^2}\right)}{\sqrt{2 \pi \sigma_i^2}}$$ and $p(z)$ is also normal. I did search the site and I think that answers to n-fold symbolic integral in Mathematica might be applicable but I really couldn't understand what was going on. Any clues are welcome, Thanks.

UPDATE: I tried to manually compute the integral for product of 3 terms by the following expression $$\int \frac{\sqrt{2 \pi }}{e^{z^2} \left(\sqrt{2 \pi \sigma _1^2} e^{\frac{\left(x_1-z\right){}^2}{2 \sigma _1^2}}\right) \left(\sqrt{2 \pi \sigma _3^2} e^{\frac{\left(x_3-z\right){}^2}{2 \sigma _3^2}}\right) \left(\sqrt{2 \pi \sigma _2^2} e^{\frac{\left(x_2-z\right){}^2}{2 \sigma _2^2}}\right)} \, dz$$

And Mathematica produced a pretty bad expression: $$\tiny\frac{\sigma _1 \sigma _2 \sigma _3 \exp \left(-\frac{\left(\left(2 \sigma _3^2+1\right) \sigma _2^2+\sigma _3^2\right) x_1^2-2 \sigma _2^2 x_3 x_1+\left(\left(2 \sigma _2^2+1\right) \sigma _1^2+\sigma _2^2\right) x_3^2-2 x_2 \left(\sigma _1^2 x_3+\sigma _3^2 x_1\right)+\left(\left(2 \sigma _3^2+1\right) \sigma _1^2+\sigma _3^2\right) x_2^2}{2 \left(\left(\left(2 \sigma _3^2+1\right) \sigma _2^2+\sigma _3^2\right) \sigma _1^2+\sigma _2^2 \sigma _3^2\right)}\right) \text{erf}\left(\frac{\sigma _1^2 \left(\sigma _2^2 \left(-x_3+2 \sigma _3^2 z+z\right)+\sigma _3^2 \left(z-x_2\right)\right)+\sigma _2^2 \sigma _3^2 \left(z-x_1\right)}{\sqrt{2} \sigma _1 \sigma _2 \sigma _3 \sqrt{\left(\left(2 \sigma _3^2+1\right) \sigma _2^2+\sigma _3^2\right) \sigma _1^2+\sigma _2^2 \sigma _3^2}}\right)}{2 \sqrt{2 \pi } \sqrt{\sigma _1^2} \sqrt{\sigma _2^2} \sqrt{\sigma _3^2} \sqrt{\left(\left(2 \sigma _3^2+1\right) \sigma _2^2+\sigma _3^2\right) \sigma _1^2+\sigma _2^2 \sigma _3^2}}$$

Answer 1

We can manual expand the symbolic sum with the function linearExpand given in my answer to How to do algebra on unevaluated integrals?.

Clear[linearExpand];
linearExpand[e_, x_, head_] := 
  e //. {op : head[arg_Plus, __] :> Distribute[op], 
    head[arg1_Times, rest__] :> 
     With[{dependencies = Internal`DependsOnQ[#, x] & /@ List @@ arg1},
      Pick[arg1, dependencies, False] head[Pick[arg1, dependencies, True], rest]
      ]};

One more difficulty to overcome is that Sum is evaluated throughout (perhaps) the integration process. It really slows things down. The expansion of the Sum results in sums that are constant -- the z can be factored out of each sum. We need to inactivate the sum. But for some reason Inactivate does not make the Sum Inactive, and wrapping the head Sum in Inactive does not speed up the integration. So let's resort to manual intervention and replace Sum by an inert head sum. Then the integration finishes in a few seconds.

exponent = linearExpand[
  MapAt[Apart[#, z] &, Sum[-(z - x[i])^2/(2 σ[i]^2), {i, n}], 1], i, Sum]

ClearAll[sum];
SetAttributes[sum, HoldAll];
With[{int = Exp[exponent /. {Sum -> sum}]},
 Integrate[int, {z, -Infinity, Infinity}, 
   Assumptions -> (Coefficient[exponent, z^2] < 0 /. Sum -> sum)] /. sum -> Sum
 ]

Answer 2

First note that since the $p(z)$ and the $p(x_i|z)$ are all Gaussians, the expression can be rewritten as $$I=\int \prod_{i=1}^n p(x_i | z) p(z) dz=A\int\prod_{i=1}^{2n}\exp\left(-c_i^2(z-z_i)^2\right)=A\int\exp\left(-\sum_{i=1}^{2n}c_i^2(z-z_i)^2\right)\\=A\exp(-d)\int\exp\left(-c_0^2(z-z_0)^2\right)\\=\frac{A\sqrt{\pi}\exp(-d)}{c_0}$$ where the $z_0,c_0,d$ arise from completing the square of the resulting quadratic polynomial $\sum_{i=1}^{2n}c_i^2(z-z_i)^2$. Mathematica can instantly factor such expressions, so this should be fairly easy to implement.

Here's an example, using $A=1,z_i=c_i=i$ and a code snippet from J.M that completes the square:

$c = Table[i, {i, 10}];
$z = Table[i, {i, 10}];
$s = Sum[$c[[i]]^2 (z - $z[[i]])^2, {i, 10}];
depress[poly_] := depress[poly, First@Variables[poly]]
depress[poly_, x_] /; PolynomialQ[poly, x] := 
 Module[{n = Exponent[poly, x], x0}, 
  x0 = -Coefficient[poly, x, n - 1]/(n Coefficient[poly, x, n]);
  Normal[Series[poly, {x, x0, n}]]]
depress[$s]

which produces

$$\sum_{i=1}^{10}c_i^2(z-z_i)^2=385 \left(z-\frac{55}{7}\right)^2+\frac{10956}{7}$$

which means that

$$I=\frac{\sqrt{\pi}\exp\left(-\frac{10956}{7}\right)}{\sqrt{385}}$$

which is an astronomically tiny number, as expected.