Computing probabilities from joint distribution of functions of random variables

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!

[bym]

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)

Answer 1

[Roger Stafford]

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.

[Steve Essinger]

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]

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])