n <- 1000
m <- n/4
x <- rnorm(n,0,1)
x1 <- x[1:m]
x2 <- x[(m+1):(2*m)]
x3 <- x[(2*m+1):(3*m)]
x4 <- x[(3*m+1):(4*m)]
x <- cbind(x1,x1+2*x2,-2*x1+3*x2-x3,x2-x3+x4)

mean(x)
sd(x)

par(mfrow=c(2,2))
qqnorm(x1)
abline(a=0,b=1)
qqnorm(x2)
abline(a=0,b=1)
qqnorm(x3)
abline(a=0,b=1)
qqnorm(x4)
abline(a=0,b=1)

xbar <- colMeans(x)
S <- var(x)
D12 <- diag(sqrt(1/diag(V)))
R <- D12%*%S%*%D12
R

S.inverse <- solve(S)
S%*%S.inverse
S.inverse%*%S

dist.M <- rep(0,m)

for (i in 1:m) {
	dist.M[i] <- t(x[i,]-xbar)%*%S.inverse%*%(x[i,]-xbar)	}
sample.quantiles <- sort(dist.M)

p <- seq(from=1/(2*m),by=1/m,length=m)
theoretical.quantiles <- qchisq(p, df=4)
par(mfrow=c(1,1))

plot(theoretical.quantiles,sample.quantiles)
abline(a=0,b=1)


MNQuantilePlot <- function(x) {

	xbar <- colMeans(x)
	S <- var(x)

	S.inverse <- solve(S)
	dist.M <- rep(0,m)

	for (i in 1:m) {
		dist.M[i] <- t(x[i,]-xbar)%*%S.inverse%*%(x[i,]-xbar)	}
	sample.quantiles <- sort(dist.M)

	p <- seq(from=1/(2*m),by=1/m,length=m)
	theoretical.quantiles <- qchisq(p, df=4)
	par(mfrow=c(1,1))

	plot(theoretical.quantiles,sample.quantiles,xlab="Theoretical quantiles",
		ylab="Sample quantiles")
	abline(a=0,b=1)
}

MNQuantilePlot(x)