Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Computing probabilities from joint distribution of functions of random variables

Asked by Steve Essinger on 24 Jan 2013

I have three independent random variables, each from a normal distribution with zero mean, but respective variances.

x1 ~ N(0,s1)

x2 ~ N(0,s2)

x3 ~ N(0,s3)

I need to find the joint probability: P(x1+x3>0, x2-x3 >0)

My approach has been to form functions of the random variables, create a joint distribution and then find the probabilities. However, I'm having difficulty coding this correctly in Matlab. Any thoughts?

My functions of random variables:

y1 = x1+x3

y2 = x2-x3

y3 = x3 % included because I need as many equations as unknowns

So,

x1 = y1-y3

x2 = y2-y3

x3 = y3

The Jacobian of the transformation is equal to 1.

So,

fY = 1/((2*pi)^(3/2)*det(SIGMA)^(1/2))*exp(-0.5*Y'*inv(SIGMA)*Y),

where Y = [y1-y3, y2-y3, y3] and SIGMA = [s1,0,0; 0,s2,0; 0,0,s3]

Now I need to compute P(y1>0, y2>0). Any ideas? Best practices? Thanks!

1 Comment

bym on 29 Jan 2013

since your means are 0, then the P(x1+x3>0) = .50 The same would hold true for x2-x3>0. Then I would say the answer of both being true is .25 (unless I am missing something entirely)

Steve Essinger

1 Answer

Answer by Roger Stafford on 24 Jan 2013

You can use matlab's 'integral3' function with numerical integration taken with respect to y1, y2, and y3. Your limits of integration would be:

 y1min = 0, y1max = inf, y2min = 0, y2max = inf, y3min = -inf, y3max = inf

Your Jacobian is 1 so you don't have to include it in the integrand. The integrand would be obtained by substituting in your joint density function using x1, x2, and x3, the corresponding y1, y2, and y2 values in accordance with the transformation you have worked out. Note however that your equation for x2 is incorrect and should be

 x2 = y2+y3

and not the one you have written.

2 Comments

Steve Essinger on 28 Jan 2013

Thank you, Roger. Good eye on the sign for equation for x2.

I'm using the integral3 function on version R2012a, but keep receiving errors on the implementation:

Error using integralCalc/finalInputChecks (line 515) Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the 'ArrayValued' option to true.

Ok, so I add this parameter to the function call:

p = integral3(fun,y1min,y1max,y2min,y2max,y3min,y3max,'ArrayValued','true');

and I receive this error:

Error using integral3 (line 109) Argument 'ArrayValued' did not match any valid parameter of the parser.

I have tried other combinations of 'ArrayValued', 'true' pairs, but nothing seems to work.

My test function is:

fun = @(y1,y2,y3) 1/((2*pi)^(3/2)*det(SIGMA)^(1/2))*exp(-0.5*[y1;y2;y3]'*inv(SIGMA)*[y1;y2;y3]);

which if I evaluate separately plugging in values for y1,y2,y3 I get a single number as expected. Note that I get the same set of errors if I change my function to something more basic:

fun = @(y1,y2,y3) [y1;y2;y3]'*[y1;y2;y3];

How do I resolve this ArrayValued issue?

bym on 29 Jan 2013

I don't have the latest Matlab, so I don't know if this will help but you could try:

vectorize(fun)
ans =
@(y1,y2,y3)1./((2.*pi).^(3./2).*det(SIGMA).^(1./2)).*exp(-0.5.*[y1;y2;y3]'.*inv(SIGMA).*[y1;y2;y3])
Roger Stafford

Contact us