Return to Question

I've been led to believe (see here and here) that Mahalanobis distance is the same as the Euclidean distance on the PCA-rotated data. In other words, taking multivariate normal data $X$, the Mahalanobis distance of all of the $x$'s from any given point (say $\mathbf{0}$) should be the same as the Euclidean distance of the entries of $X^{rot}$ from $\mathbf{0}$, where $X^{rot}$ is the product of the data and the PCA rotation matrix.

I've been led to believe (see here and here) that Mahalanobis distance is the same as the Euclidean distance on the PCA-rotated data. In other words, taking multivariate normal data $X$, the Mahalanobis distance of all of the $x$'s from any given point (say $\mathbf{0}$) should be the same as the Euclidean distance of the entries of $X^{rot}$ from $\mathbf{0}$, where $X^{rot}$ is the product of the data and the PCA rotation matrix.

Two questions:

Two questions:

f = function(N=1000,n,p){
    a = runif(p^2,-1,1)
    a = matrix(a,p)
    S = t(a)%*%a
    x = mvrnorm(N,rep(0,p),S)
    mx = apply(x, 2, mean)
    sx = apply(x, 2, sd)
    x = t(apply(x,1,function(X){(X-mx)/sx}))
    Ss = solve(cov(x))
    x = x[sample(1:N,n,replace=F),]
    md = mahalanobis(x,rep(0,p),Ss,inverted=T)
    prMD<-pchisq(md,df = p) 
    pc = prcomp(x,center=F,scale=F)
    d<-mahalanobis(scale(pc$x),rep(0,ncol(pc$x)),diag(rep(1,ncol(pc$x))))
	prPCA<-pchisq(d,df = min(p,n))#N is the number of PCs where N<P 
	return(data.frame(prbias = as.numeric(mean(prMD - prPCA)), prvariance = as.numeric(mean((prMD - prPCA)^2))))
}
grid = data.frame(n=100,p=2:200)
grid$prvariance <-grid$prbias <-NA
for (i in 1:nrow(grid)){
	o = f(n=grid[i,]$n,p=grid[i,]$p)
    grid[i,3:4]<-o
}
par(mfrow=c(1,2))
with(grid, plot(p,prbias))
abline(v=100)
m = lm(prbias~p,data=grid)
abline(m,col='red',lty=2)
with(grid, plot(p,prvariance))
abline(v=100)

f = function(N=1000,n,p){
  a = runif(p^2,-1,1)
  a = matrix(a,p)
  S = t(a)%*%a
  x = mvrnorm(N,rep(0,p),S)
  mx = apply(x, 2, mean)
  sx = apply(x, 2, sd)
  x = t(apply(x,1,function(X){(X-mx)/sx}))
  Ss = solve(cov(x))
  x = x[sample(1:N,n,replace=F),]
  md = mahalanobis(x,rep(0,p),Ss,inverted=T)
  prMD<-pchisq(md,df = p) 
  pc = prcomp(x,center=F,scale=F)
  d<-mahalanobis(scale(pc$x),rep(0,ncol(pc$x)),diag(rep(1,ncol(pc$x))))
  prPCA<-pchisq(d,df = min(p,n))#N is the number of PCs where N<P 
  return(data.frame(prbias = as.numeric(mean(prMD - prPCA)), prvariance = as.numeric(mean((prMD - prPCA)^2))))
}

grid = data.frame(n=100,p=2:200)
grid$prvariance <-grid$prbias <-NA

for (i in 1:nrow(grid)){
    o = f(n=grid[i,]$n,p=grid[i,]$p)
    grid[i,3:4]<-o
} 

par(mfrow=c(1,2))
with(grid, plot(p,prbias))
abline(v=100)
m = lm(prbias~p,data=grid)
abline(m,col='red',lty=2)
with(grid, plot(p,prvariance))
abline(v=100)