Convolution between two PDFs using product of Laplace transforms in Symbolic Toolbox
6 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
Paul
le 26 Fév 2023
Modifié(e) : Paul
le 27 Fév 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 commentaires
Voir également
Catégories
En savoir plus sur Assumptions dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!