## R-code (September 20, 2007)

## univariate normal commands

# draw 100 realizations from a normal random variable with mean 5, and variance 4.
mu <- 5
std.dev <- sqrt(4)
y.sam <- rnorm(100, mean=mu, sd=std.dev)

summary(y.sam)

par(mfrow=c(2,1))
plot(y.sam)
plot(density(y.sam), lwd=2)

# the density function is dnorm()
# let's draw 100 samples and plot the true density
mu <- 0
std.dev <- sqrt(1)
y.sam <- rnorm(100, mean=mu, sd=std.dev)

x <- seq(-4,4, by=0.1)
f.x <- dnorm(x, mean=mu, sd=std.dev)
f.y.sam <- dnorm(y.sam,mean=mu, sd=std.dev)

plot(y.sam, f.y.sam, pch=16, col="blue", xlim=c(-4,4))
lines(x, f.x, type="l", lwd=2)
lines(density(y.sam), col="blue", lwd=2)


## for a normal random variable with mean 0 and std.dev=1, what is:

## P(y < 0.5)
pnorm(0.5, 0, 1)

## p(y > 0.5)
1-pnorm(0.5, 0, 1)

## let's see this
x <- seq(-4,4, by=0.1)
f.x <- dnorm(x, 0, 1)
plot(x, f.x, type="l", lwd=2)
abline(v=0.5, lwd=2, col="red")

## multivariate normal commands
## make two plots with samples from a draw from a population with 1.) a low value for the det of Sigma 2.) higher value.
Mu <- matrix(c(0,0),2,1)
Sigma.low <- matrix(c(1,0.9,0.9,1),2,2)
Sigma.high.1 <- matrix(c(1,0.3,0.3,1),2,2)
Sigma.high.2 <- matrix(c(10,0.9,0.9,10),2,2)

det(Sigma.low)
det(Sigma.high.1)
det(Sigma.high.2)

## let's look at the eigen values
eigen(Sigma.low)$values
eigen(Sigma.high.1)$values
eigen(Sigma.high.2)$values

## load MASS library for mvrnorm()
library(MASS)

par(mfrow=c(2,2))
y <- mvrnorm(100, Mu, Sigma.low)
plot(y[,1], y[,2], pch=16)

y <- mvrnorm(100, Mu, Sigma.high.1)
plot(y[,1], y[,2], pch=16)

y <- mvrnorm(100, Mu, Sigma.high.2)
plot(y[,1], y[,2], pch=16)


## in class assignment.  write a function called my.mvdnrom(y, Mu, Sigma) to calculate the density of a multivariate random variable:
## here is a start:

my.mvdnorm <- function(y, Mu, Sigma){
out <- ## fill in density
return(out)
}

## there is a library called mvtnorm, which has this function so le't test it.
Sigma <- matrix(c(5,2,2,3),2,2)
Mu <- matrix(c(2,4), 2,1)
y <- matrix(c(0,0), 2,1)

## install and download package from web and load library for mvtnorm()
library(mvtnorm)
dmvnorm(t(y), Mu, Sigma)

my.mvdnorm(y, Mu, Sigma)