Clear Filters
Clear Filters

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

7 views (last 30 days)
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.

Answers (1)

Paul
Paul on 26 Feb 2023
Edited: Paul on 27 Feb 2023
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).

Products


Release

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!