The $p$ variate normal distribution is parameterized by a $p$ dimensional mean vector $\mu$ and a $p \times p$ variance-covariance matrix $\Sigma$.

We say ${\bf Y} =(Y_1, \ldots, Y_p)^T$ is a random vector variate that follows a multivariate normal distrbution, ${\bf Y} \sim N_p(\mu, \Sigma)$, and has density: $$f({\bf Y}) = f(Y_1, \ldots, Y_p) = \left(2 \pi \right)^{-\frac{p}{2}} det(\Sigma)^{-{\frac{1}{2}}} e^{-\frac{1}{2} \left({\bf Y} -\mu \right)^T \Sigma^{-1} \left({\bf Y} -\mu \right)}$$

When $p=1$ this is the usual normal density and we write $Y \sim N(\mu, \sigma^2)$.

It might be helpful to consider what this density looks like in particular situations. We’ll use the package mvtnorm in R.

require('mvtnorm')

## Loading required package: mvtnorm

We will plot contours from a few multivariate normals when $p=2$.

#
# First get a grid 
#
xmin <- -10
xmax <- 10
ymin <- -10
ymax <- 10

x <- seq(xmin, xmax, by = 0.1)
nxvals <- length(x)
y <- seq(ymin, ymax, by = 0.1)
nyvals <- length(y)
gridLocs <- cbind(rep(x, times=nyvals), rep(y, each= nxvals))
indexes <- cbind(rep(1:nxvals, times=nyvals), rep(1:nyvals, each= nxvals))
density <- matrix(0, nrow=nxvals, ncol=nyvals)

# Construct a common (empty) plot
# 
newplot <- function() {
  plot(0, type="n", xlim=c(xmin,xmax),
     ylim=c(ymin,ymax), 
     asp=1, 
     xlab="x", ylab="y" )
}

We will plot data from a few multivariate normals when $p=2$.

# 

newplot()
# First standard multivariate normal contours
mu_std <- c(0,0)
sigma_std  <- diag(1, 2, 2)
density[indexes] <-  dmvnorm(gridLocs, mean = mu_std , sigma=sigma_std )
contour(x=x, y=y, z=density, add=TRUE)

# relocate
density[indexes] <-  dmvnorm(gridLocs, mean = c(5,5) , sigma=sigma_std )
contour(x=x, y=y, z=density, col="red", add=TRUE)

# relocate and double standard deviation of x
density[indexes] <-  dmvnorm(gridLocs, mean = c(-5,5) , 
                             sigma=matrix(c(4,0,
                                            0,1),
                                          nrow=2, byrow=TRUE))
contour(x=x, y=y, z=density, col="blue", add=TRUE)

# relocate and double standard deviation of y
density[indexes] <-  dmvnorm(gridLocs, mean = c(-5,-5) , 
                             sigma=matrix(c(1,0,
                                            0,4),
                                          nrow=2, byrow=TRUE))
contour(x=x, y=y, z=density, col="green", add=TRUE)

# relocate, double standard deviation of one, and rotate by 45 degrees = pi/4
theta <- pi/4
rot <- matrix(c( cos(theta), -sin(theta),
                 sin(theta),  cos(theta)),
              nrow=2, byrow=TRUE)
sigma <- matrix(c(1,0,
                        0,4),
                      nrow=2, byrow=TRUE)

sigma <- t(rot) %*% sigma %*% rot
density[indexes] <-  dmvnorm(gridLocs, mean = c(5,-5),
                             sigma = sigma
                             )
contour(x=x, y=y, z=density, col="purple", add=TRUE)

Note you could also generate points from these distributions

# A common plot
newplot()
# First standard bivariate normal
y <- rmvnorm(500, mean=mu_std, sigma_std )
points(y, col="black", pch=19,cex=0.5)


# relocate
y <- rmvnorm(500, mean = c(5,5) , sigma=sigma_std )
points(y, col="red", pch=19,cex=0.5)

# relocate and double standard deviation of x
y <- rmvnorm(500,  mean = c(-5,5) , 
             sigma=matrix(c(4,0,
                            0,1),
                          nrow=2, byrow=TRUE))
points(y, col="blue", pch=19,cex=0.5)

# relocate and double standard deviation of y
y <- rmvnorm(500, mean = c(-5,-5) , 
             sigma=matrix(c(1,0,
                            0,4),
                          nrow=2, byrow=TRUE))
points(y, col="green", pch=19,cex=0.5)

# relocate, double standard deviation of one, and rotate by 45 degrees = pi/4
theta <- pi/4
rot <- matrix(c( cos(theta), -sin(theta),
                 sin(theta),  cos(theta)),
              nrow=2, byrow=TRUE)
sigma <- matrix(c(1,0,
                  0,4),
                nrow=2, byrow=TRUE)

sigma <- t(rot) %*% sigma %*% rot
y <- rmvnorm(500, mean = c(5,-5), sigma = sigma)
points(y, col="purple", pch=19,cex=0.5)