CODEWITHSUNDEEP
matlabparametersnumerical-integration
Alex
Alex

Reputation: 603

Need Help in Writting Matlab Code for Sum of Independent Random Variables

I want to write a MATLAB code for the convolution of density functions .

For random variables X (with PDF f_x)and Y(with PDF f_y), the PDF of their summation, X+Y, can be obtained by the following MATLAB code:

function p=sumX_Y(z,mean1,sigma1,mean2,sigma2)
G = @(y)f_x(mean1,sigma1,z-y).*f_y(mean2,sigma2,y);p= integral(G,-Inf,Inf);

Now for summation of three random variable X+Y+Z it should be 

function p=sumX_Y_Z(z,mean1,sigma1,mean2,sigma2,mean3,sigma3)
    G = @(y)sumX_Y(z-y,mean1,sigma1,mean2,sigma2).*f_z(mean3,sigma3,y);p= integral(G,-Inf,Inf);

But it will not work since the first input of the function sumX_Y is not of type double. I would appreciate if you could help me to fix this problem.

I also want to find the PDF for sum of more than 3 random variable and I do not know how to write the algorithm. Thanks a lot!

Upvotes: 1

Views: 1247

Answers (2)

Buck Thorn
Buck Thorn

Reputation: 5073

One way to perform the convolution is as follows:

First define sumX_Y

function p=sumX_Y(z,mean1,sigma1,mean2,sigma2)

f_x = @(x) normpdf(x,mean1,sigma1);
f_y = @(x) normpdf(x,mean2,sigma2);


for ii=1:length(z)

    G = @(y)f_x(z(ii)-y).*f_y(y);
    p(ii)= quad8(G, ...
                 min([mean1,mean2])-3*max([sigma1,sigma2])  ,...
                 max([mean1,mean2])+3*max([sigma1,sigma2])    );
end

The loop in the function above can be vectorized using the bsxfun notation btw.

Then define sumX_Y_Z as follows:

function p=sumX_Y_Z(z,mean1,sigma1,mean2,sigma2,mean3,sigma3)

f_z = @(x) normpdf(x,mean3,sigma3);
G = @(y)sumX_Y(z-y,mean1,sigma1,mean2,sigma2).*f_z(y);

p= quad8(G,min([mean1,mean2,mean3])-3*max([sigma1,sigma2,sigma3]),...
           max([mean1,mean2,mean3])+3*max([sigma1,sigma2,sigma3]));

The whole thing can be called as follows:

mean1=1;
sigma1=0.3;
mean2=2;
sigma2=1;
mean3=0;
sigma3=1;

X = 1.4; % <-- value at which to evaluate p

p=sumX_Y(X,mean1,sigma1,mean2,sigma2)
p=sumX_Y_Z(X,mean1,sigma1,mean2,sigma2,mean3,sigma3)

with output:

p =

    0.1181


p =

    0.1494

You can modify the functions to use integral instead of 'quad8', and substitute your definition of f_z etc of course.

Upvotes: 2

devrobf
devrobf

Reputation: 7213

You can convert a number to a double like this:

numberAsDouble = double(number);

So either convert the number before calling sumX_Y, or force it to be converted in the function itself:

function p=sumX_Y(z,mean1,sigma1,mean2,sigma2)
  G = @(y)f_x(mean1,sigma1,double(z)-y).*f_y(mean2,sigma2,y);p= integral(G,-Inf,Inf);

You could do this for all the inputs if you like.

Upvotes: 0

Related Questions