# Plotting the convolution of the same probability density function with itself

I would like to compute the integral

$$q_n(x)=\int_{B}^A q_{n-1} (\omega)f(x-\omega)\mathrm{d} \omega,\quad q_1=f,\quad n\geq 1.$$

for a given density function $f$. Basically it is the convolution of a pdf with itself for $n$ times in the range $A$ and $B$. when $n=1$, I simply plot the original pdf $f$, if $n=2$, then I convolve it with itself only once and plot etc..

Eventually I would like to plot the densities $q_n$ for $n=1,2,3..,7$ in the same figure.

My code with the help of some people from this forum is

f0[x_] := PDF[NormalDistribution[-2, 2], x]
f1[x_] := PDF[NormalDistribution[2, 2], x]

q[1, B_, A_, f_] := f

q[n_ /; n > 1, B_, A_, f_][x_] := Module[{\[Omega]},
tempIntegrate[Evaluate[q[n - 1, B,A, f]][\[Omega]] f[x - \[Omega]], {\[Omega],B, A}]]

p[n_, B_, A_, f_] :=
tempIntegrate[q[n, B, A, f][x], {x, -\[Infinity], B}]+tempIntegrate[q[n, B, A, f][x], {x,A, \[Infinity]}] //.
{s_ tempIntegrate[b_, a__] :>tempIntegrate[s b, a],
tempIntegrate[tempIntegrate[b_, a__], c__] :>tempIntegrate[b, a, c]} /.
tempIntegrate -> NIntegrate


I can not

Plot[q[2, -2, 2, f0][x],{x,-10,10}]


correctly and I cannot see the analytical formulation when I type

q[2, B, A, f0][x]


EDIT :

My eventual intention was to calculate

$$P[n]=\int_{(-\infty,B)\cup (A,\infty)} q_n (x) \mathrm{d}x$$

and p[2, -2, 2, f0] provides me a numerical value but I am also interested in some analytical values.

To do this I tried to change ->NIntegrate to Integrate but then the code didnt function properly.

• Maybe you are not finished editing yet, but what is your definition of tempIntegrate? – Jacob Akkerboom Oct 10 '13 at 12:36
• If you replace tempIntegrate with Integrate, then Plot[q[2, 1, 3, f0][y], {y, 0, 2}] works for me. Note that Evaluate does nothing in your code. Evaluate is often used in Plot to ignore the fact that Plot does not evaluate its first argument. Integrate however does not hold any arguments. In addition, the Evaluate would only override this if it would surround the entire first argument. I suggest removing the Evaluate :). Note that in the Plot code I said that worked for me, the function gets integrated a lot of times. I think that is not necessary – Jacob Akkerboom Oct 10 '13 at 12:43
• @JacobAkkerboom I tried your suggestions before but the results which I got was incorrect. Please give me 5 minutes and I get some correct numerical values for comparison. – Seyhmus Güngören Oct 10 '13 at 12:50
• 0.5227501319481026 , 0.2973909442488435 , 0.11261520561472473 , 0.04211024363961929 , 0.01573950408880347 are the true values for $p[n, -2, 2, f0]$ for $n=1,...,5$ – Seyhmus Güngören Oct 10 '13 at 12:55
• @JacobAkkerboom i did your suggestion, removed evaluate and changed tempintegrate to Integrate and eventually p[3, -2, 2, f0] gave me nothing)) – Seyhmus Güngören Oct 10 '13 at 13:08

Maybe this will do what you want

ClearAll[q]
q[1, B_, A_, f_] := f
q[n_ /; n > 1, B_, A_, f_] :=
Block[{\[Omega]},
q[n, B, A, f] =
Function[x,
Evaluate[
Integrate[
q[n - 1, B, A, f][\[Omega]] f[x - \[Omega]], {\[Omega], B, A}]]
]
]


Plotting q2 will now be very fast, as Mathematica will evaluate the integral for q2 to a simpler form, even before x is known, which greatly speeds up computation. That is, in this case you do not evaluate an integral for every x you wish to plot.

However, it turns out Mathematica is not able to find a nice expression for q3. We can still plot it, but at every step Mathematica tries to evaluate a difficult integral analytically, so this takes long and is suboptimal. Maybe I can improve it :).

Example:

Plot[q[2, 1, 3, f0][x], {x, 0, 2}]


Note that the definition of q[2,1,3,f0] is now

q[2,1,3,f0]

(*output*)
Function[x$,(E^(-(1/16) (4+x$)^2) (-Erf[(2-x$)/4]+Erf[(6-x$)/4]))/(8 Sqrt[\[Pi]])]


so that we indeed see that here the Integral has been analytically "solved"/"simplified"/"evaluated".

However, again, this does not work for q3. The definition of that becomes

q[3, 1, 3, f0]
(*output in InputForm ;)*)

Function[x$, Integrate[(E^(-(2 + x$ - \[Omega])^2/8 - (4 + \[Omega])^2/16)*
(-Erf[(2 - \[Omega])/4] + Erf[(6 - \[Omega])/4]))/(16*Sqrt[2]*Pi), {\[Omega], 1, 3}]]


where we see that the function Integrate is still present.

• Note that to get values for q4 you have to use the function N, that is, write for example q[4, 1, 3, f0][2] // N. – Jacob Akkerboom Oct 10 '13 at 13:21

There may be a better way to approach this problem. If you transform your function f[t] into the frequency domain ft[w], convolution becomes multiplication, and the iterative convolution of f[t] with itself becomes raising ft[w] to a power. Here is a simple example where

f[t_] := Exp[-t^2/2];
ft[w_] := FourierTransform[f[t], t, w]


Now you can plot the powers of ft[w]:

Plot[Table[ft[w]^n, {n, 1, 5, 1}], {w, -3, 3}]


which are the transforms of the time functions. You can return to the time domain by inverse-transforming, for example, at the 5th step,

InverseFourierTransform[ft[w]^5, w, t]

E^(-(t^2/10))/Sqrt[5]


Of course you can replace this f[t] with your f0[t] without any problems. Below, the OP asked for how to do this when the range of the function is restricted. This can be handled by truncating the function. Here's a truncated normal distribution:

f1[t_, a_, b_] := UnitStep[t - a] UnitStep[b - t] PDF[NormalDistribution[-2, 2], t];


This can be transformed (almost) as easily as above:

f1t[w_, a_, b_] = FourierTransform[f1[t, a, b], t, w];
Plot[Table[N[Abs[f1t[w, -3, 3]]]^n, {n, 1, 5, 1}], {w, -3, 3}, PlotRange -> All]


• Thanks for the answer. It is not a bad idea but my convolutions is in some prespecified range $A$ and $B$. how can I do in in a parametric way? – Seyhmus Güngören Oct 10 '13 at 15:00
• Define your function to be whatever you want in the region A to B, and zero outside. UnitStep is one way to do this. – bill s Oct 10 '13 at 15:01
• thank you very much for the edit. I have bigger problem but the idea of frequency domain sounds very interesting, After some progress I will get back to you. – Seyhmus Güngören Oct 10 '13 at 15:43
• I dont get the result of inverse Fourier transform – Seyhmus Güngören Oct 15 '13 at 20:21
con[f_, g_, x_, l_, u_] :=
Integrate[(f /. x -> x ) (g /. x -> (y - x)), {x, l, u}];
func[u_, ll_, ul_] :=
Nest[con[#, PDF[NormalDistribution[], x], x, ll, ul] &,
PDF[NormalDistribution[], x], u];


Using integration limits of -1 and 1:

tab = Table[func[j, -1, 1], {j, 5}];
Plot[tab, {y, -3, 3}, PlotLegends -> Range[5]]


gives:

EDIT

To allow use of f0[x] and f1[x] and incorporating alternative suggested by bills):

q[f_, g_, x_, b_, a_, n_] :=
Nest[Convolve[#, f (UnitStep[x - b] - UnitStep[x - a]), x, y] &, g,
n]


Applying to f0[x] (note must put as expression of x not just f0]:

qtab = Table[q[f0[x], f0[x], x, -1, 1, j], {j, 5}];


and plotting:

Plot[qtab, {y, -4, 4}, PlotLegends -> Range[5]]


• It might be easier to use the Convolve command, which does more or less what your con does. – bill s Oct 11 '13 at 13:01
• @bills thank you for the feedback. Convolve uses -Infinity to Infinity as limits of integration (as far as I understand the documentation). I could not find options to alter the limits of integration, hence custom function. – ubpdqn Oct 11 '13 at 13:14
• You can limit the range of integration by changing the function (Unitstep` for instance) because any given convolution has fixed endpoints. But your custom function is fine too... – bill s Oct 11 '13 at 14:16
• @bills thank you for this, another option for my future reference – ubpdqn Oct 11 '13 at 14:17
• I am not able to use it in my program. I replaced PDF[NormalDistribution[] parts with either $f_0$ or $f_0[x]$ and both gave an error. I would like to have a form q[n,B,A,f] for which I can select an $f$ too, for example $f_0$ as given in the question. How one should modify the code? – Seyhmus Güngören Oct 11 '13 at 14:46