Bivariate normal value standardization

I want to standardize a bivariate normal CDF. I tried with inverse square root of covariance matrix and with Cholesky decomposition. The results are always different across all 3. I don't know why.

sigma=[1,0.5;0.5,1];
X = [1,1];
z=mvncdf(X,[0,0],sigma);

%%method 1
X1=X*sqrtm(inv(sigma));
z1=mvncdf(X1,[0,0],[1,0;0,1]);
%method 2
L = chol(sigma, 'lower');
X11=X*inv(L);
z11=mvncdf(X11,[0,0],[1,0;0,1]);
%results

disp([z,z1,z11])

    0.7452    0.6287    0.5814

[Umar]

Hi CJ,

To ensure consistency in standardization, you can try using a standardized input vector by transforming X using the mean and standard deviation of the bivariate normal distribution.

Answer 1

[Paul]

Hi CJ,

In short, the area of integration for the X1 case is no longer a rectangle as is assumed by mvncdf

Define the original distribution of an MVN vector X

Sigma = [1,0.5;0.5,1];
mu = [0 0];

Find the probability that -inf < X1 < 1 & -inf < X2 < 1

X = [1,1];
p1 = mvncdf(X,mu,Sigma)

p1 = 0.7452

This probablity can also be computed by integrating under the pdf of X.

Find the pdf of X

x = -3:.01:3;
[X1,X2] = meshgrid(x);
pdf1 = reshape(mvnpdf([X1(:),X2(:)],mu,Sigma),size(X1));

Plot it and add the limits at X1 = 1 and X2 = 1;

figure
pcolor(X1,X2,pdf1),shading interp
xline(1,'w');yline(1,'w');

We can approximate the probability by numerical integration under the pdf over the lower left square of the plot. Of course we are not capturing the tails of the density.

mask = (X1 <= 1) & (X2 <= 1);
trapz(x,trapz(x,pdf1.*mask,2),1)

ans = 0.7443

Same (close enough) result as above.

Let Z be standard MVN, we have X = A*Z, where

A = sqrtm(Sigma)

A = 2x2

0.9659 0.2588 0.2588 0.9659

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Plot the pdf of Z

z = -3:.01:3;
[Z1,Z2] = meshgrid(z);
pdf2 = reshape(mvnpdf([Z1(:),Z2(:)],mu,eye(2)),size(Z1));
figure
pcolor(Z1,Z2,pdf2),shading interp,colorbar

Now, to properly compute the probability we need to find the region in the Z-plane that maps through

X = A*Z

to the lower left square above in the X-plane.

% X = A*Z
mask = reshape(all((A*[Z1(:),Z2(:)].').' <= [1 1],2),size(Z1));
hold on

Overlay the mask on the Z-plane to visualize the region of integration (which extends down and left to infinity)

plotmask = double(mask);
plotmask(plotmask == 0) = nan;
scatter3(Z1(1:10:end,1:10:end),Z2(1:10:end,1:10:end),plotmask(1:10:end,1:10:end),'w.'),view(2)

Compute the probability in z-space

trapz(z,trapz(z,pdf2.*mask,2),1)

ans = 0.7426

[CJ]

That makes sense, thank you very much! It doesn't seem like there's a way to use MVNCDF then to do what I want? I have to calculate the integral manually? I was hoping to do this transformation on a vector many times and mvncdf is quite efficient (i.e. provide many X and Sigmas, transform, evaluate on standard normal as vector). Any suggestions?

[Paul]

You're very welcome.

As far as I know, mvncdf can only be used over rectangular regions, possibly extending to -inf in two directions.

I'm not sure what the issue is. If you have a non-standard normal vector, like X above, and want to find the probability over a rectangular region, why not just use mvncdf? Why transform to a standard normal vector?

[CJ]

I have to evaluate the CDF millions of times, which is time consuming. I have code that closely approximates an uncorrelated normal (from the bvnl function) that is 3x faster than mvncdf. Hence I want to transform a correlated into an uncorrelated one.