# Transform general multivariate normal to standard multivariate normal

I have data that I can assume will be multivariate normal with a known mean vector mu and known covariance matrix sigma, and I'm looking to identify points that fall outside the 95% ellipse. I think the proper way to do this is to apply a transform to the data to a standard multivariate normal, compute the euclidean distances to zero, and compare against the value derived from running mvtnorm::qmvnorm(.95,tail='both',mean=mu,sigma=sigma). However, if this is correct (and feel free to correct me if there's a simpler solution to labelling the points outside the 95% ellipse), I'm struggling with the transform part. I know I'd subtract the mean vector mu from the observed data matrix, and I need to do something involving the resulting matrix and sigma, but I'm lost on what. Some code below showing some data and ending where I'm stuck:

mu = c(100,50)
sigma = matrix(c(100,210,210,900),nrow=2)
obs_data = MASS::mvrnorm(1e2,mu=mu,Sigma=sigma)
transformed = t(t(obs_data) - mu)
#what now?

• – whuber
Aug 19, 2020 at 18:44
• Ah! So the Mahalanobis distance is the euclidean distance-from-zero after the transform, so I can just use the stats::mahalanobis() function. Neat! Quick follow-up @whuber: do compare these distances to the 95% ellipse (or hypersphere) would I be getting the comparison value via mvtnorm::qmvnorm(.95,tail='both',mean=mu,sigma=sigma) or mvtnorm::qmvnorm(.95,tail='both',sigma=diag(2))? (i.e. do I need to tell qmvrnorm the original mean & covariance?) I'm getting very different results between these two. Aug 19, 2020 at 20:03
• Oh, looking at the help page for stats::mahalanobis() and it shows that the distances should follow the same quantile function as a chi-square distribution with df=ncol(x), so I don't even need qmvnorm. When I use qchisq(.95,df=ncol(dat)), I do indeed observe that 95% of the distances exceed this value. So I wonder what qmvnorm() is doing or is supposed to be used for? Aug 19, 2020 at 20:57
• qmvnorm gives you the half-width of a mean-centered square (or, generally, hypercube) in the original coordinates. In the transformed coordinates that region would correspond to a generalized rhombus.
– whuber
Aug 19, 2020 at 21:02

With the help of the content here, I was able to work this out and thought I'd post the resulting code and verification:

#define a function to transform a standard multivariate normal (zero mean,
# identity covariance) to a non-standard multivariate normal with specified
# means and covariance.
std_to_mvn = function(Z,mu,sigma){
tchol_sigma = t(chol(sigma))
X_0centered = t(tchol_sigma %*% t(Z))
return(X)
}

#define a function to transform a multivariate normal with mean mu and covariance
# sigma to a standard multivariate normal with zero mean and identity covariance
mvn_to_std = function(X,mu,sigma){
tchol_sigma = t(chol(sigma))
X_0centered = t(t(X)-mu)
Z = t(solve(tchol_sigma,t(X_0centered)))
}

#Check that we can induce the non-standard multivariate normal

#first,sample from a "standard" multivariate normal
z = MASS::mvrnorm(
n = 1e5
, mu = rep(0,times=3)
, Sigma = diag(3)
)

#next, specify the parameters (arbitrary; taken from https://jamesmccaffrey.wordpress.com/2017/11/03/example-of-calculating-a-covariance-matrix/)
mu = c(68, 600, 40)
sigma = matrix(
c(
11.5 , 50, 34.75
, 50, 1250, 205
, 34.75, 205, 110
)
, nrow = 3
, ncol = 3
)

#now, transform and check that we get what we expected
x = std_to_mvn(z,mu,sigma)
max(abs(colMeans(x)-mu)) # <1
max(abs(cov(x)-sigma)) # <1

#Check that we can transform a non-standard multivariate normal to standard
x = MASS::mvrnorm(
n = 1e5
, mu = mu
, Sigma = sigma
)

#now, transform and check that we get what we expected
z = mvn_to_std(x,mu,sigma)
max(abs(colMeans(z))) # <1
max(abs(cov(z)-diag(3))) # <1