Convolution between two PDFs using product of Laplace transforms in Symbolic Toolbox
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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.
0 Kommentare
Antworten (1)
Paul
am 26 Feb. 2023
Bearbeitet: Paul
am 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
F(x) = 2.*(1 - 0.01)*dirac(x) + rectangularPulse(0, 400, x).*0.01./400
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)
But Matlab can't seem to evaluate that integral for some reason. But we can do this
int(E(x),x,-1e10,1e10)
Or this
expand(int(E(x),x,-inf,inf))
And we see that E(x) is not a valid pdf. Same for F(x).
0 Kommentare
Siehe auch
Kategorien
Mehr zu Assumptions finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!