# Convolution of a denstity function having a point mass with itself

My density function is as follows:

der[x_] :=Piecewise[{{E^((-4 - 4*x - x^2 - Log^2)/8)/(2^(x/4)*Sqrt[Pi]*(2*Sqrt[2/E]*Erf[Log/(2*Sqrt)]+Erfc[(-2 + Log)/(2*Sqrt)] +2*Erfc[(2 + Log)/(2*Sqrt)])), x > 0}, {0.2015152271993863, x == 0}},(2^(x/4)*E^((-4 - 4*x - x^2 - Log^2)/8))/(Sqrt[Pi]*(2*Sqrt[2/E]*Erf[Log/(2*Sqrt)]+Erfc[(-2 + Log)/(2*Sqrt)] + 2*Erfc[(2 +Log)/(2*Sqrt)]))]


This density has a point mass of

0.2015152271993863


at $x=0$.

Taking into account the point mass, I would like to calculate $P[n];$

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

with

$$q_n(x)=\int_{B}^A q_{n-1} (\omega)f(x-\omega)\mathrm{d} \omega,\quad q_1=f,\quad n\geq 1.$$ and $f$ is some density function,which is $\mathrm{der}[x]$ here. To make these iterative calculations, I have the following code:

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


If the calculations are all correct, then the summation of all elements of the following table should (almost) add up to $1$

tab2 = Table[p[j, -2, 2, der], {j, 10}];


Since I cannot deal with the point mass it doesnt add upto $1$. If you want to see that it indeed adds up to $1$, you can consider for example

der[x_]:=PDF[NormalDistribution[-2, 2], x]


How can I deal with the point mass in the convolution and also in the integrations?

## 1 Answer

One way to deal with the point mass is to use the function DiracDelta to represent the discontinuity. Using your density function der[x], the integral is:

c = Integrate[der[x], {x, -Infinity, Infinity}]

0.798485


So the height of the delta function must be 1-c (=0.201515). Accordingly, your PDF is:

pdf[x_]:=der[x] + (1 - c) DiracDelta[x]


To check, we can integrate and make sure the integral is unity:

Integrate[pdf[x], {x, -Infinity, Infinity}]


which is indeed 1.

• yes the density is becoming a real density but my iterations are still not working. – Seyhmus Güngören Oct 17 '13 at 14:14
• This doesnt help me to solve my problem, dont you have any comments about the convolution? – Seyhmus Güngören Oct 17 '13 at 16:00
• I suggested here: mathematica.stackexchange.com/a/33799/1783 how you might solve this, but you chose to ignore my suggestion. Using Fourier transforms is a standard trick for bypassing and simplifying convolutions. Above, I have answered your question about how to include the point mass in the pdf. – bill s Oct 17 '13 at 16:22
• I didnt ignore. I run all the answers given by all the persons. I did it one by one. Fourier idea is a good one and it can really be helpful to my whole problem. That is the reason why I left it to the end. I run exactly your suggestion and the values that I got for p[n_, B_, A_, f_] were incorrect. Here are the true results: 0.5227501319481026 , 0.2973909442488435 , 0.11261520561472473 , 0.04211024363961929 , 0.01573950408880347 for p[n,−2,2,f0] for n=1,...,5. With your suggestion I didnt get the true results. I think the problem is about the "range of integration", $(B,A)$ – Seyhmus Güngören Oct 17 '13 at 21:49
• the range $(B,A)$ applies to the result of the convolution $f_0*f_0$. This means we apply the range to the output of the convolution not to the individual densities. As a result I got incorrect results. You can also check but I think very likely that the results that you will obtain with the current codes that you suggested will not be correct either. I checked your other solution which you gave as a comment to the other answer too(about the command 'convolve'). It gives the same incorrect result. For $n=3$, instead of the correct result 0.11261520561472473` I got 0.13... – Seyhmus Güngören Oct 17 '13 at 22:05