Using convolution to determine pdf of adding two triangular random variables

Hello,

In my research, I'm representing a likelihood of investment success with triangular distribution.

I want to know the likelihood of two investments succeeding given amount y, or P(x1+x2 <y), where x1 and x2 are different triangular distributions. So far my search has shown me that I need to use conv function, which I did below:

clear all;
clc;

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

x = linspace(85,210,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,'valid');

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;

There are two problems:

1. PDF of X1+X2 integrates to much higher than 1, given by variable p12.
2. PDF of X1+X2 varies widely depending on the range of X I use.

Intuitively I feel that the valid range of X1+X2 should begin at 85, since probability of investment succeeding is zero for either x1 or x2 at this point. Also, shouldn't X1+X2 end at 110+100, the sum of upper limit of x1 and x2?

So what is wrong with my code? Help would be greatly appreciated!

Answer 1

[Robert]

The conv function is essentially doing a numerical integration over x. It is not necessary that both pdfs be based on the same x, but both x arrays must have the same constant spacing between values. I think you will find that the following code does what you want:

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

    dx = 0.01; % to ensure constant spacing
    x1 = 85:dx:100; % Could include some of the region where
    x2 = 90:dx:110; % 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);

    pdf1 = pdf(pd1,x1);
    pdf2 = pdf(pd2,x2);

    pdf12 = conv(pdf1,pdf2)*dx;

    p1  = trapz(x1,pdf1)
    p2  = trapz(x2,pdf2)
    p12 = trapz(x12,pdf12)

    plot(x1,pdf1,x2,pdf2,x12,pdf12)