Using integral3 to calculate the conditional expectation of an RV

Hi,

I am trying to calculate E[X|X>Y & X>Z], where X,Y,Z~N(0,1). I am trying to use the function integral3, however I am unable to get the correct answer. Here is my code:

xmin = -inf;
xmax = inf;
ymin = -inf;
ymax = @(x) x;
zmin = -inf;
zmax = @(x,y) x;
integral3(@(x,y,z) x.*normpdf(x).*normpdf(y).*normpdf(z),xmin,xmax,ymin,ymax,zmin,zmax)

The answer I get, 0.282, is exactly a third of what the correct answer is: 0.846, which I calculate the following way:

eps = normrnd(0,1,10000000,3);
temp = eps(:,1)>eps(:,2) & eps(:,1)>eps(:,3);
mean(eps(temp,1),1)

Could someone advise me what I'm doing wrong?

P.S. The reason I don't simply use the latter method is that, if I generalize the condition to E[X|X>Y+a & X>Z+b], where a,b are constants, then for large a,b occurrences when X>Y+a and X>Z+b become extremely rare, and the method breaks down

Answer 1

[Roger Stafford]

Since you are computing the conditional expected value, you need to divide your integral by the probability of satisfying X>Y & X>Z. That is, divide by

integral3(@(x,y,z) normpdf(x).*normpdf(y).*normpdf(z),xmin,xmax,ymin,ymax,zmin,zmax)

Your approximation by "mean(eps(temp,1),1)" does that automatically in the 'mean' operation.

Note that you have assumed the independence of X, Y, and Z without explicitly stating it.

Also note that approximating with 10000000 random numbers will give you much less than seven decimal-place accuracy in your result. Performing the numerical integration is by far the most accurate.

[Krzysztof]

You're absolutely right, thank you!