I am wonder how to standardize multivariate normal value.
Normal standard multivariate distribution of $q$ variables is $z\sim N_q(0,I_q)$. I have found that $a+Bz\sim N_q(Ba,BB^T)$ and based on this fact normalization could be performed throught two ways:
1) Substracting mean vector and then taking $B=\Sigma^{-0.5}$ (as it gives $\Sigma^{-0.5} \Sigma (\Sigma^{-0.5})^T=I_q$) we get normalized value: $z=\Sigma^{-0.5}(x-a)$.
2) Using Cholesky decomposition $\Sigma=U^TU$ we get normalization rule $z=(U^T)^{-1}(x-a)$
But both strategies does not work!
So you can check it using matlab code below:
MU=[0.1333,-0.1];
SIGMA=[0.95,0.19;0.19,0.97];
x=[0.5, 0.7];
mvncdf(x,MU,SIGMA)
z=(SIGMA^(-0.5))*transpose(x-MU);
mvncdf(z)
z=transpose(chol(SIGMA))^(-1)*transpose(x-MU);
mvncdf(z)
I also add code for R that gives the same results proving that the problem is not due to programm but due to incorrect appoach:
SIGMA=matrix(c(0.95,0.19,0.19,0.97), ncol=2)
MU=c(0.1333, -0.1)
x=c(0.5,0.7)
pmvnorm(mean = MU, sigma=SIGMA, lower=-Inf, upper=x)
z=as.numeric((solve(sqrtm(SIGMA)))%*%(x-MU))
pmvnorm(mean = MU, sigma=SIGMA, lower=-Inf, upper=z)
As the difference is rather big I am sure that such standardization approach is incorrect. But it is the only technicks I have found despite it is a hot and applied subject.
Will be very greatful for help!
PS: I also performed Jordan form approach but the result is still wrong:
[v,j]=jordan(SIGMA)
B=v*j^(-0.5)*v^(-1)
z=B*transpose(x-MU);
mvncdf(z)
Best Answer
If $X\sim N_q(a,\Sigma)$ (so $X$ is a random vector taking values in $\mathbb R^{q\times1}$) and $B\in\mathbb R^{p\times q}$ then $BX\sim N_p(Ba,B\Sigma B^T)$.
The matrix $\Sigma^{-1/2}$ should be taken to be the positive-definite symmetric square root of $\Sigma^{-1}$. Then we have $\Sigma^{-1/2}(X-a)\sim N_q(0,I_q)$.
Where I'm sitting I don't have access to MATLAB. Maybe if you could post some output I could say something about that. My initial suspicion is that MATLAB may be taking $\Sigma^{-1/2}$ to be something other than the positive-definite symmetric square root of $\Sigma^{-1}$. If you could tell us what MATLAB tells you is the value of $\Sigma^{-1/2}$ that would tell us whether that is what is happening.