CODEWITHSUNDEEP
matlabstatisticsprobabilityconvolutionprobability-distribution
Yohan K.
Yohan K.

Reputation: 23

calculating sum of two triangular random variables (Matlab)

I would like to calculate the sum of two triangular random variables,

P(x1+x2 < y)

image

Is there a faster way to implement the sum of two triangular random variables in Matlab?

EDIT: It seems there's possibly a much easier way, as shown in this minitab demonstration. So it's not impossible. It doesn't explain how the PDF was calculated, sadly. Still looking into how I can do this in matlab.

EDIT2: Following advice, I'm using conv function in Matlab to develop the PDF of the sum of two random variables:

clear all;
clc;


pd1 = makedist('Triangular','a',85,'b',90,'c',100);
pd2 = makedist('Triangular','a',90,'b',100,'c',110);

x = linspace(85,290,200);
x1 = linspace(85,100,200);
x2 = linspace(90,110,200);
pdf1 = pdf(pd1,x1);
pdf2 = pdf(pd2,x2);

z = median(diff(x))*conv(pdf1,pdf2,'same');

p1  = trapz(x1,pdf1) %probability P(x1<y)
p2  = trapz(x2,pdf2) %probability P(x2<y)
p12 = trapz(x,z)     %probability P(x1+x2 <y)

hold on;
plot(x1,pdf1) %plot pdf of dist. x1
plot(x2,pdf2) %plot pdf of dist. x2
plot(x,z)     %plot pdf of x1+x2
hold off;

However this code has two problems:

  1. PDF of X1+X2 integrates to much higher than 1.
  2. PDF of X1+X2 varies widely depending on the range of x. Intuitively, if the X1+X2 is larger than 210 (the sum of upper limits "c" of two individual triangular distributions, 100 + 110), shouldn't P(X1+X2 <210) equal to 1? Also, since the lower limits "a" is 85 and 90, P(X1+X2 <85) = 0?

Upvotes: 1

Views: 1520

Answers (2)

Yohan K.
Yohan K.

Reputation: 23

Below is the proper implementation for future users. Many thanks to Robert Dodier for guidance.

clear all;
clc;

min1 = 85;
max1 = 100;
min2 = 90;
max2 = 110;
y    = 210;

pd1 = makedist('Triangular','a',min1,'b',90,'c',max1);
pd2 = makedist('Triangular','a',min2,'b',100,'c',max2);

dx = 0.01; % to ensure constant spacing
x1 = min1:dx:max1; % Could include some of the region where
x2 = min2:dx:max2; % the pdf is 0, but we don't have to.

x12 = linspace(...
    x1(1)   + x2(1) , ...
    x1(end) + x2(end) , ...
    length(x1)+length(x2)-1);
[c,index] = min(abs(x12-y));
x_short = linspace(min1+min2,x12(index),index);

pdf1 = pdf(pd1,x1);
pdf2 = pdf(pd2,x2);
pdf12 = conv(pdf1,pdf2)*dx;

zz = pdf12(1:index);
zz(index) = 0;

p1  = trapz(x1,pdf1)
p2  = trapz(x2,pdf2)
p12 = trapz(x_short,zz)

plot(x1,pdf1,x2,pdf2,x12,pdf12)
hold on;
fill(x_short,zz,'blue')     % plot x1+x2
hold off;

Upvotes: 1

Robert Dodier
Robert Dodier

Reputation: 17575

The pdf of the sum of independent variables is the convolution of the pdf's of the variables. So you need to compute the convolution of two variables with trianular pdf's. A triangle is piecewise linear, so the convolution will be piecewise quadratic.

There are a few ways to about it. If a numerical result is acceptable: discretize the pdf's and compute the convolution of the discretized pdf's. I believe there is a function conv in Matlab for that. If not, you can take the fast Fourier transform (via fft), compute the product point by point, then take the inverse transform (ifft if I remember correctly) since fft(convolution(f, g)) = fft(f) fft(g). You will need to be careful about normalization if you use either conv or fft.

If you must have an exact result, the convolution is just an integral, and if you're careful with the limits of integration, you can figure it out by hand. I don't know if the Matlab symbolic toolbox is available to you, and if so, I don't know if it can handle integrals of functions defined piecewise.

Upvotes: 2

Related Questions