I am trying to perform the symbolic integration of the following expression

∫n∏i=1p(xi|z)p(z)dz\int \prod_{i=1}^n p(x_i | z) p(z) dz

where p(xi|z)=exp(−(x−z)22σ2i)√2πσ2ip(x_i|z) = \frac{\exp\left(\tfrac{-(x-z)^2}{2\sigma_i^2}\right)}{\sqrt{2 \pi \sigma_i^2}} and p(z)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

∫√2πez2(√2πσ21e(x1−z)22σ21)(√2πσ23e(x3−z)22σ23)(√2πσ22e(x2−z)22σ22)dz\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:

σ1σ2σ3exp(−((2σ23+1)σ22+σ23)x21−2σ22x3x1+((2σ22+1)σ21+σ22)x23−2×2(σ21×3+σ23×1)+((2σ23+1)σ21+σ23)x222(((2σ23+1)σ22+σ23)σ21+σ22σ23))erf(σ21(σ22(−x3+2σ23z+z)+σ23(z−x2))+σ22σ23(z−x1)√2σ1σ2σ3√((2σ23+1)σ22+σ23)σ21+σ22σ23)2√2π√σ21√σ22√σ23√((2σ23+1)σ22+σ23)σ21+σ22σ23\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}}

=================

If I’m understanding this correctly, note that ∏ni=1p(xi|z)p(z)dz\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.

– DumpsterDoofus

Nov 13 ’14 at 1:39

Also, could you explain what the multiple-integral tag is for? Your expression appears to be the integral over a single variable zz.

– DumpsterDoofus

Nov 13 ’14 at 1:40

@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.

– Pushpendre

Nov 13 ’14 at 1:43

Removed the multiple integral tag

– Pushpendre

Nov 13 ’14 at 1:44

=================

2 Answers

2

=================

First note that since the p(z)p(z) and the p(xi|z)p(x_i|z) are all Gaussians, the expression can be rewritten as

I=∫n∏i=1p(xi|z)p(z)dz=A∫2n∏i=1exp(−c2i(z−zi)2)=A∫exp(−2n∑i=1c2i(z−zi)2)=Aexp(−d)∫exp(−c20(z−z0)2)=A√πexp(−d)c0I=\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 z0,c0,dz_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\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=iA=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}\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}}I=\frac{\sqrt{\pi}\exp\left(-\frac{10956}{7}\right)}{\sqrt{385}}

which is an astronomically tiny number, as expected.

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.

– Pushpendre

Nov 13 ’14 at 2:15

@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 nn is a bit hard to do in Mathematica. Could you give an example of what you would want to do?

– DumpsterDoofus

Nov 13 ’14 at 2:18

We can manual expand the symbolic sum with the function linearExpand given in my answer to

How to do algebra on unsolved 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

]