1

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?
SecretAgentMan
  • 2,856
  • 7
  • 21
  • 41
Yohan K.
  • 23
  • 5
  • To perform an element wise multiplication you need a `.*` between the two pdfs. So `fun = @(x) pdf(pd2,x).*pdf(pd1,y-x);` – schvaba986 Dec 21 '15 at 09:35
  • Regarding terminology, the probability function you define depends on a single argument, so I wouldn't call it a 'joint' probability. – Itamar Katz Dec 21 '15 at 10:30
  • Edited for incorrect understanding on subject. Sorry! The question is now properly phrased. – Yohan K. Dec 21 '15 at 15:36

2 Answers2

2

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.

Robert Dodier
  • 16,905
  • 2
  • 31
  • 48
  • Thank you for your guidance. It seems `conv` is the right way to go about this in Matlab. I've updated my question to show the progress so far. – Yohan K. Dec 22 '15 at 08:52
  • Looks like you have got most of the solution. You have to be careful about the normalization of the pdf's in order to make sure that the pdf's of x, y, and x + y all integrate to 1. Also, note that the support (nonzero density) of x + y is the interval (min(x) + min(y), max(x) + max(y)) which you have to map onto the result from `conv`. Specifically, the first element of the `conv` result corresponds to the the first element of x + first element of y, likewise last element corresponds to last of x + last of y. There are length(x) + length(y) - 1 elements in the `conv` result. – Robert Dodier Dec 22 '15 at 19:11
1

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;
Yohan K.
  • 23
  • 5