1
$\begingroup$

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}}$$

$\endgroup$
4
  • $\begingroup$ If I'm understanding this correctly, note that $\prod_{i=1}^n p(x_i | z) p(z) dz$ is also a Gaussian after rearrangement of the terms in the exponent, so you can bypass any such symbolic difficulties by first factoring the integrand, and then doing the single integral of a single Gaussian, which Mathematica can do easily. I'll post an answer in a while. $\endgroup$ Nov 13, 2014 at 1:39
  • $\begingroup$ Also, could you explain what the multiple-integral tag is for? Your expression appears to be the integral over a single variable $z$. $\endgroup$ Nov 13, 2014 at 1:40
  • $\begingroup$ @DD: Actually I had a more complicated integral in mind but simplified it to the guassian example. I'll remove the multiple integral tag since its not applicable to the simple case that I posted. $\endgroup$
    – Pushpendre
    Nov 13, 2014 at 1:43
  • $\begingroup$ Removed the multiple integral tag $\endgroup$
    – Pushpendre
    Nov 13, 2014 at 1:44

2 Answers 2

4
$\begingroup$

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
 ]

Mathematica graphics

$\endgroup$
3
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Thanks for the answer DD, just one one more clarification though. Could you mention how to symbolically represent a summation of n terms when n is not known beforehand ? I am a mathematica noob and I couldn't find this in the documentation. $\endgroup$
    – Pushpendre
    Nov 13, 2014 at 2:15
  • $\begingroup$ @Pushpendre: I actually don't know how to do that, unfortunately; Mathematica is usually used for manipulating concrete expressions where the terms are provided. Manipulating symbolic sums for arbitrary $n$ is a bit hard to do in Mathematica. Could you give an example of what you would want to do? $\endgroup$ Nov 13, 2014 at 2:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.