#################################################################
#### Illustration of Multivariate Gaussian ######################
#### Anup Aprem #################################################
#################################################################
library(fields)

# Mean of multivariate Gaussian
mu  	=  	c(1, 0) 

# Covariance matrix f
S  	=  	matrix(c(4, 3,3, 4), nrow=2, ncol=2, byrow=TRUE) # covariance

# Eigen Decomposition of Covariance matrix
E  	=  	eigen(S)
D  	=  	diag(E$values)
U  	=  	E$vectors

# Generate grid to plot multivariate Gaussian
stepSize = 	 0.5
x  	=  	seq(-5,5,stepSize)
y  	=  	seq(-5,5,stepSize)
m 	= 	length(x); n=length(y)
X 	= 	matrix(rep(x,each=n),nrow=n)
Y 	= 	matrix(rep(y,m),nrow=n)
r  	=  	dim(X)[1]
l  	=  	dim(X)[2]
data  	=  	cbind(as.vector(X),as.vector(Y))

# Evaluate the multivariate density
a  	=  	(2*pi)^(-(2/2)) * det(S)^(-0.5)
b  	=  	matrix(0,1,dim(data)[1])
for (ij in 1:dim(data)[1]) {
	u  	=  	(data-mu)[ij,]
	b[ij]   =  	exp(-.5 * t(u) %*% solve(S) %*% (u))
}
p= a * b;
dim(p) = c(r, l)
p=p*(stepSize^2)

# Plot the density and the countour lines
attach(mtcars)
par(mfrow=c(2,1))
persp(predictSurface(Tps(as.matrix(expand.grid(x=x,y=y)),as.vector(p))),main='Plot of multivariate Gaussian',zlab = "Density",xlab='x',ylab='y')
contour(predictSurface(Tps(as.matrix(expand.grid(x=x,y=y)),as.vector(p))),nlevels=10)
title(main='Contour plot of multivariate Gaussian',xlab='x',ylab='y')
# Plot first eigenvector
lines(c(mu[1], mu[1]+sqrt(D[1,1])*U[1,1]),c(mu[2], mu[2]+sqrt(D[1,1])*U[2,1]),lwd=3,col='red')
# Plot second eigenvector
lines(c(mu[1], mu[1]+sqrt(D[2,2])*U[1,2]),c(mu[2], mu[2]+sqrt(D[2,2])*U[2,2]),lwd=3,col='red')