Example R programs and commands
13. Estimating, sampling and plotting bivariate densities
# All lines preceded by the "#" character are my comments.
# All other left-justified lines are my input.
# All other indented lines are the R program output.
##### ESTIMATING the mean and covariance of a bivariate normal from samples
# Read 17 samples from a bivariate normal density:
# x1 x2 i (of 17 samples)
#---------- ---------- --
data <- scan()
1.033024 0.418857 1
0.322703 0.613913 2
-1.949918 0.505329 3
-3.910608 0.081451 4
-1.163915 0.904262 5
-4.944068 -0.086801 6
0.936824 0.889408 7
-2.697000 0.185941 8
-1.484696 0.318464 9
-2.969133 -0.019170 10
-0.565761 0.682739 11
0.314755 0.010014 12
3.181993 0.862641 13
-0.816990 0.440558 14
-2.053569 0.444235 15
-0.311087 0.797796 16
-1.131430 0.522798 17
# Impose structure: 3 columns (x1, x2, i) with 17 rows of data triplets.
x123 <- matrix(data, ncol=3, byrow=T); x123
# NOTE: R fills columns first by default, in the top-to-bottom then
# left-to-right order. Use byrow=TRUE (or byrow=T) to fill by rows.
# Extract the first two (data) columns, omitting the third (index i) column,
# giving default names to the rows and columns:
x<- data.frame(x123[,1:2]); x
# Calculate the mean of the columns:
eMu <- colMeans(x); eMu
# Compute the sum-of-squares matrix and divide by the number of degrees
# of freedom (number of samples less one) to get the covariance matrix:
# using R's covariance function:
eSig <- cov(x); eSig
# Find the eigenvalues (and eigenvectors) of the covariance matrix:
eigen(eSig)
## The result is estimates for the mean (mu) and covariance (Sigma)
## parameters of the bivariate normal density that produced the
## original samples.
##### GENERATE random samples with the estimated mean and covariance:
## Use the function mvrnorm() after adding the MASS package to R:
library("MASS") # perhaps after 'install.packages("MASS")'
n <- 17 # since there are 17 old samples
sim <- mvrnorm(n, mu=eMu, Sigma=eSig);
all <- rbind(x,sim) # concatenate the old and new points.
blackred <- gl(2,n,2*n) # color labels: 1=black, 2=red
plot(all, col=blackred) # compare old and new samplings
## Alternative: obtain random samples with the same covariance directly by
## factoring the covariance matrix to find the right rotation:
library("mgcv") # load the matrix square root function from its library
y<- rbind(rnorm(n), rnorm(n)) # two rows of indep N(0,1) coordinates
new <- mroot(eSig)%*%y + eMu # mix the coordinates to get new points
new <- t(new) # ...and transpose to get one new point per row
all <- rbind(x,new) # concatenate the old and new points.
blackred <- gl(2,n,2*n) # color labels: 1=black, 2=red
plot(all, col=blackred) # compare old and new samplings
# The covariance of the 'new' samples should approximate 'cov(x)' but
# is almost certain to be different:
cov(new) # covariance of the new points
cov(x) # covariance of the old points
##### PLOT a bivariate normal density from its formula.
## Use "contour()" and "persp()" after producing samples
# Use the function "bvnpdf()" defined on the class website under LINKS
source("bvnpdf.R") # after downloading into your R folder,...
#... or cut and paste it into your R session
# Generate x,y coordinates for a (square) grid of points to evaluate:
x<-seq(-9,9,by=0.25); y<-x;
# Evaluate "bvnpdf()" at the grid of points:
z<-bvnpdf(x,y,mu=eMu, Sigma=eSig );
# Perspective plot of the graph in (x,y,z):
persp(x,y,z)
# Contour plot of the graph in (x,y,z):
contour(x,y,z)