Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

My density function is as follows:

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

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?

share|improve this question
add comment

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.

share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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