# Multivariate normal value standardization

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)


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.

• Matlab gives me $\Sigma^{-0.5}$=(0.7207,-0.0587;-0.0587,1.0235). Initial matrix is positive defined as egivenvalues are 0.9345 and 1.9855. Real value I anticipate is 0.4934 while after standardization I get 0.4614 and for other examples the difference could be even greater. Commented Feb 13, 2016 at 0:10
• I have also checked that $\Sigma^{-1/2}$ is positive defined Commented Feb 13, 2016 at 0:21
• Working in R, if I let $B$ be the matrix $\left[\begin{array}{rr} 0.7207 & -0.0587; \\ -0.0587 & 1.0235 \end{array} \right]$, then the command solve(B)%*%solve(B), which gives the value of $(B^{-1})^2$, yields $\left[\begin{array}{rr} 1.9497729 & 0.1899401 \\ 0.1899401 & 0.9699799 \end{array} \right]$ rather than $\left[\begin{array}{rr} 0.95 & 0.19 \\ 0.19 & 0.97 \end{array} \right]$. So something is wrong with your way of trying to find $\Sigma^{-1/2}$. $\qquad$ Commented Feb 13, 2016 at 0:25
• For $\Sigma^{-1/2}$, I'm getting $\left[ \begin{array}{rr} 1.0415497 & -0.1035578 \\ -0.1035578 & 1.0306489 \end{array} \right]$. $\qquad$ Commented Feb 13, 2016 at 0:29
• Sorry I have mistype. I had 1.95 instead of 0.95 but it is do not change the problem substantialy. Finding $\Sigma^{-1/2}$ in my approach is correct as because multiplying it then both sides on $\Sigma$ gives me identical vector. Commented Feb 13, 2016 at 0:36

Before I came to this site, using Wolfram Mathematica to calibrate a gyroscope (MPU6050), I came up with the same result as you simply intuiting that multi-dimensional spread might be worth. $\text{gRdist}\text{:=}\text{EstimatedDistribution}\left[\text{data},\text{MultinormalDistribution}\left[\left\{\mu _1,\mu _2,\mu _3\right\},\left( \begin{array}{ccc} v_{11} & v_{12} & v_{13} \\ v_{21} & v_{22} & v_{23} \\ v_{31} & v_{32} & v_{33} \\ \end{array} \right)\right]\right];$

$gYoffset = Mean[gRdist];$

$gYscale = Inverse[MatrixPower[Covariance[gRdist], 1/2]];$

$\{-42.0838,149.375,28.0686\}$

$\left( \begin{array}{ccc} 0.0833029 & 0.000290934 & \text{3.567851038398982$\grave{ }$*${}^{\wedge}$-7} \\ 0.000290934 & 0.0745029 & 0.000244209 \\ \text{3.567851038399714$\grave{ }$*${}^{\wedge}$-7} & 0.000244209 & 0.0818224 \\ \end{array} \right)$

$\text{adjustedData}=(\text{gYscale}.(\text{$\#$1}-\text{gYoffset})\&)\text{/@}\text{data}$

$\text{EstimatedDistribution}\left[\text{adjustedData},\text{MultinormalDistribution}\left[\left\{\mu _1,\mu _2,\mu _3\right\},\left( \begin{array}{ccc} v_{11} & v_{12} & v_{13} \\ v_{21} & v_{22} & v_{23} \\ v_{31} & v_{32} & v_{33} \\ \end{array} \right)\right]\right]$

WONDERFULLY!!!

$\text{MultinormalDistribution}\left[\{\text{2.0686693247474613$\grave{ }$*${}^{\wedge}$-16},\text{7.8355722154314$\grave{ }$*${}^{\wedge}$-16},\text{1.7542859740454608$\grave{ }$*${}^{\wedge}$-16}\},\left( \begin{array}{ccc} 1. & -\text{3.309187133831902$\grave{ }$*${}^{\wedge}$-16} & \text{1.9671686402361793$\grave{ }$*${}^{\wedge}$-17} \\ -\text{3.2239708270537383$\grave{ }$*${}^{\wedge}$-16} & 1. & \text{3.7404711706817533$\grave{ }$*${}^{\wedge}$-16} \\ \text{1.9671686402361793$\grave{ }$*${}^{\wedge}$-17} & \text{3.7404711706817533$\grave{ }$*${}^{\wedge}$-16} & 1. \\ \end{array} \right)\right]$

If I use //Chop, the result is:

$\text{MultinormalDistribution}\left[\{0,0,0\},\left( \begin{array}{ccc} 1. & 0 & 0 \\ 0 & 1. & 0 \\ 0 & 0 & 1. \\ \end{array} \right)\right]$

Thanks!