Convolution between two PDFs using product of Laplace transforms in Symbolic Toolbox

Hello.

I am trying to obtain probability density function (PDF) H(x) resulting from the the convolution between two probability density functions E(x) and F(x) - the objective is to obtain the PDF of the sum between two random variables. I want to do it using the Symbolic Toolbox, calculating the convolution as the product of the Laplace transforms of E(x) and F(x). Here is what I am doing.

>> syms x

>> E(x) = 2.*(1 - 0.01)*dirac(x) + rectangularPulse(0, 300, x).*0.01./300;

>> F(x) = 2.*(1 - 0.01)*dirac(x) + rectangularPulse(0, 400, x).*0.01./400;

I am checking that the integral of these functions is 1, as should be for a PDF (integrating up to 1000 is enough):

>> int(E,x,0,1000)

ans = 1

>> int(F,x,0,1000)

ans = 1

If I calculate the Laplace transform of one of these PDFs and then apply the inverse Laplace transform, I get the same function, as expected:

>> ilaplace(laplace(E,x,s),s,x)

ans = (99*dirac(x))/50 - heaviside(x - 300)/30000 + 1/30000

>> int(ilaplace(laplace(E,x,s),s,x), 0, 1000)

ans = 1

OK, now I generate the Laplace transform of the convolution between E(x) and F(x), as G(s):

>> syms s

>> LE=laplace(E,x,s)

LE = 1/(30000*s) - exp(-300*s)/(30000*s) + 99/50

>> LF=laplace(F,x,s)

LF = 1/(40000*s) - exp(-400*s)/(40000*s) + 99/50

>> G=LE*LF

G = (1/(30000*s) - exp(-300*s)/(30000*s) + 99/50)*(1/(40000*s) - exp(-400*s)/(40000*s) + 99/50)

Now, I generate the resulting PDF, by inverting the Laplace transform:

>> H=ilaplace(G,s,x)

H = x/1200000000 - (33*heaviside(x - 300))/500000 - (99*heaviside(x - 400))/2000000 + (9801*dirac(x))/2500 - (heaviside(x - 300)*(x - 300))/1200000000 - (heaviside(x - 400)*(x - 400))/1200000000 + (heaviside(x - 700)*(x - 700))/1200000000 + 231/2000000

And this is where I find the problem. H(x) is not a PDF, since its integral is greater than 1:

>> int(H, x, 0, 1000)

ans = 19999/10000

Is there someone able to tell me what I am doing wrong and how I should correct it? Thank you in advance.

Answer 1

[Paul]

E(x) and F(x) are not valid pdfs.

syms x real
E(x) = 2.*(1 - 0.01)*dirac(x) + rectangularPulse(0, 300, x).*0.01./300 

E(x) = 

F(x) = 2.*(1 - 0.01)*dirac(x) + rectangularPulse(0, 400, x).*0.01./400

F(x) = 

Because each has an impulse at the origin and is zero for x < 0, the integration has to start to the left of x = 0. Actually, the formal check would be that this integral is unity

int(E(x),x,-inf,inf)

ans = 

But Matlab can't seem to evaluate that integral for some reason. But we can do this

int(E(x),x,-1e10,1e10)

ans = 

Or this

expand(int(E(x),x,-inf,inf))

ans = 

And we see that E(x) is not a valid pdf. Same for F(x).