convolution of two probability density functions
显示 更早的评论
I am considering two exponential probability distribution functions with mean equal to 5 and 3.
pdf1=@(x)exppdf(x,5);
pdf2=@(x)exppdf(x,3);
I want to compute the pdf of the sum of these two densities using convolution:
x=0:100;
uh=conv(pdf1(x),pdf2(x));
Since the step-size of x is one, uh should be a density and the area under uh should be one. But
trapz(0:200,uh)
yields 1.26.
But when I choose x=0:0.1:100;
trapz(0:0.1:200,uh)*0.1
yields 1.026.
So, it turns out that the accuracy of 'using 'conv' to the get the density of the sum of two independent random variables' depends heavily upon the support chosen. Am I doing it correct or there is something wrong with my approach. Is there an automatic way to select a reasonable support for the convolution of two densities.
采纳的回答
David Goodmanson
2017-9-8
编辑：David Goodmanson
2017-9-10
Hello Abhinav
pdfs are continuous functions, so the closer spaced your x points are, the closer you get to the expected answer. In the following code the smaller you make the array spacing delx, the closer you get to the expected answer (although for delx = .0001 it takes awhile).
n1 = 3; n2 = 5;
delx = .001;
xmax = -ceil(log(1e-10)*n1)
x = 0:delx:xmax;
a1 = exp(-x/n1)/n1;
a2 = exp(-x/n2)/n2;
trapz(x,a1) % norm
trapz(x,a2.*x) % mean
trapz(x,a2)
trapz(x,a2.*x)
uh = conv(a1,a2)*delx;
x3 = 0:delx:2*xmax;
plot(x3,uh)
trapz(x3,uh)
trapz(x3,x3.*uh)
I don't have the exppdf function so I used the equivalent.
Trapz, primitive as it is, shows the trend, but it would be much preferable if Mathworks supported something better in basic Matlab like integration by cubic spline.
APPENDED
[1] Convolution of exponential pdfs
In the case of the convolution of n exponential pdfs, all of whose decay constants are unique (no repeats): let 'a' be the vector of decay constants. The kth pdf is
(1/a(k))*exp(-x/a(k)).
The convolution of all n pdfs is the sum over k of
c(k)*(1/a(k))*exp(-x/a(k)), where
c(k) = a(k)^(n-1) / [ (a(k)-a(1))*(a(k)-a(2)) ...(a(k)-a(n)) ]
In the denominator the factor (a(k) - a(k)) = 0 is excluded, so there are n-1 factors in all.
[2] Convolution by fft
For a convolution of n general pdfs, let x be a row array of N equally spaced points with spacing delx, where the range of x is wide enough that all pdfs die down to very small values at the upper and lower ends of the x array. Let M be an (n x N) matrix of the pdfs stacked on top of each other. Then the convolution is
y = real(ifft(prod(fft(M,[],2))))*delx^(n-1);
In other words take the fft of each pdf down the rows, multiply them all together down the columns and transform back. For n = 10 and N = 2e6 points this takes less than 2 seconds on my PC and it is basically linear in N (as long as N has lots of divisors).
7 个评论
Abhinav
2017-9-9
David Goodmanson
2017-9-9
Hi Abhinav,
Is there a reason you can't use the exact expression, which is
uh = (exp(-x/n1)-exp(-x/n2)) / (n1-n2) n1 ~= n2
uh = x.*exp(-x/n1) / n1^2 n1 = n2
where the pdf is 0 for x<0 just as with the exponential pdf.
Abhinav
2017-9-9
David Goodmanson
2017-9-10
Hi Abhinov,
I just modified the answer above to provide an expression for convolution of more than two exponential pdfs, and to do convolution by fft.
Abhinav
2017-9-11
Michael Reshko
2019-4-5
If we use your code
y = real(ifft(prod(fft(M,[],2))))*delx^(n-1)
to convolute, say, 10 standard normal pdfs sampled on an interval [-3, 3] , then your y does not look like a normal N(0,sqrt(10)). This code probably needs changing to take into account differences between circular and linear convolution. Are you able to provide a full solution for this?
David Goodmanson
2019-6-29
Hi Michael
I didn't attempt to reproduce what is going wrong, but there seem to be two possibilies: [1] not letting the gaussian run out far enough into the tails. You need a significant part of the array to be basically zero so that circular convolution gives the same result as linear convolution. [2] if the gaussian in x is too wide, then in the spacial frequency domain it will be very narrow and might not have enough points so that taking it to the tenth power preserves 'gaussiness'. The code below does 50 gaussians on 1000 points. I made no attempt to get the normalization totally right but you can see that the width is working correctly.
N = 1000;
x = -N/2:N/2-1;
delx = x(2)-x(1);
y = (1/sqrt(N))*exp(-4*x.^2/N);
figure(1)
plot(x,y); grid on
n = 50;
yy = ifftshift(y);
Y = repmat(yy,n,1);
ynew = real(ifft(prod(fft(Y,[],2))))*delx^(n-1);
figure(3)
plot(x,fftshift(ynew)); grid on
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
