Bivariate normal distribution with R

Exercise: Plot a bivariate normal distribution for simulated data, \mu_{x1} = 0, \sigma_{x1} = 0.5, \mu_{x2} = 0.5, \sigma_{x2}= 2, and \rho= 0.5.

As explained on Mathworld, the bivariate normal distribution is the statistical distribution with probability density function:

where:

and, with ‘real’ data:

Here, we will impose \rho= 0.5. We are free to do so as these are simulated data.

R and the S-language seem a good option for this exercise. In our code we first create two vectors of data (along the x and y dimensions). Then, then compute the densities and store them in third vector (z-dimension). We use the persp function of R to create the 3d-plot. And, that’s it.

# Édouard Tallent @ TaGoMa.Tech
# September 2012
# This code plots simulated bivariate normal distributions

# Some variable definitions
mu1 <- 0	# expected value of x
mu2 <- 0.5	# expected value of y
sig1 <- 0.5	# variance of x
sig2 <- 2	# variance of y
rho <- 0.5	# corr(x, y)

# Some additional variables for x-axis and y-axis 
xm <- -3
xp <- 3
ym <- -3
yp <- 3

x <- seq(xm, xp, length= as.integer((xp + abs(xm)) * 10))  # vector series x
y <- seq(ym, yp, length= as.integer((yp + abs(ym)) * 10))  # vector series y

# Core function
bivariate <- function(x,y){
	term1 <- 1 / (2 * pi * sig1 * sig2 * sqrt(1 - rho^2))
	term2 <- (x - mu1)^2 / sig1^2
	term3 <- -(2 * rho * (x - mu1)*(y - mu2))/(sig1 * sig2)
	term4 <- (y - mu2)^2 / sig2^2
	z <- term2 + term3 + term4
	term5 <- term1 * exp((-z / (2 *(1 - rho^2))))
	return (term5)
}

# Computes the density values
z <- outer(x,y,bivariate)

# Plot
persp(x, y, z, main = "Bivariate Normal Distribution",
		sub = bquote(bold(mu[1])==.(mu1)~", "~sigma[1]==.(sig1)~", "~mu[2]==.(mu2)~
		", "~sigma[2]==.(sig2)~", "~rho==.(rho)),
		col="orchid2", theta = 55, phi = 30, r = 40, d = 0.1, expand = 0.5,
		ltheta = 90, lphi = 180, shade = 0.4, ticktype = "detailed", nticks=5)

Notice that this code is written such that the text of the legend automatically changes accordingly with the values passed on to the  variables mu1, sig1, mu2, sig2, and rho (see at the beginning of the snippet, under # Some variable definitions). This makes things clearer and easier if the user changes the parameters to plot other simulated bivariate. Changing the values of the axes of the graph is also easy.

* exp(- z * (2 * (1 - rho^2)
About these ads

4 thoughts on “Bivariate normal distribution with R

    • Hola Alejandra.
      Gracias por visitar mi bitácora. Averigüé la linea return de mi código. Me parece correcto. Pero puede ser que ahora mis ojitos no vean nada bien porque demasiado cansancio. :)
      Saludos,
      Édouard

    • Hola Alejandra. ¡Gracías por insistir!
      Trata del term5. Hay que ahora no sirve para nada.
      Creó dicha variable pero no es imprescindible. Uno puede llevarse el resultado esperado usando return([la expresión del term5]). Non obstante la dejaré.
      Saludos,
      Édouard

Leave a Reply

Please log in using one of these methods to post your comment:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s