################################################################
# Run the following programs and print out all of the plots.
# Comment on the independence of the joint densities.

library(mvtnorm)

################################################################
# Tossing coins.

# How do the results change if the cutoff is increased to 545?
# Decreased to 515?

cutoff <- 525

# fair coins

p <- 0.5

n <- 1000

k <- 1000

x <- rbinom(k,n,p)

x.type1.err.prob <- length(x[x >= cutoff])/k  # x[x>=525] values >= 525 in x
x.type1.err.prob

# biased coins

p <- 0.55

n <- 1000

k <- 1000

x <- rbinom(k,n,p)

x.type2.err.prob <- length(x[x < cutoff])/k  # x[x<525] values < 525 in x
x.type2.err.prob

################################################################
# Sum of independent Binomial random variables.
# generate 1000 random binomials for x and y.

k <- 1000 

n <- 5
m <- 10

p <- 0.5

x <- rbinom(k, n,p)
y <- rbinom(k, m,p)

z <- x+y

X11()
hist(z,nclass=n+m)

################################################################
# Sum of independent Exponential random variables
# generate 5 random gamma random variables 1000 times

lam <- 1

k <- 1000

n <- 5

x <- matrix(rexp(k*n,lam),ncol=n)

z <- apply(x,1,sum)  # apply function to rows 1

X11()
hist(z)
X11() 
plot(density(z, n=50, window="g"),type="l",xlab="x",ylab="density")

################################################################
# Joint ditribution of U and V.  Are U and V independent?

alpha <- 5
beta <- 17
lam <- 1

k <- 1000

x <- rgamma(k, shape = alpha, rate = lam)
y <- rgamma(k, shape = beta, rate = lam)

X11()
plot(x,y)

u <- x + y 
v <- x/(x+y)

# joint

X11()
plot(u,v)

cor(u,v)

# marginals

X11()
plot(density(u, n=50, window="g"),type="l",xlab="u",ylab="density")
X11()
plot(density(v, n=50, window="g"),type="l",xlab="v",ylab="density")

################################################################
# In each case are Y_1 = X_1 + X_2 and Y_2 = X_1 - X_2 independent?
# generate 5000 random variates for x1 and x2

n <- 5000

# Uniform

x1 <- runif(n)
x2 <- runif(n)
X11()
plot(x1,x2)

y1 <- x1 + x2
y2 <- x1 - x2

X11()
plot(y1,y2)

# Exponential

lam1 <- 5
lam2 <- 10

x1 <- rexp(n, lam1)
x2 <- rexp(n, lam2)

X11()
plot(x1,x2)

y1 <- x1 + x2
y2 <- x1 - x2

X11()
plot(y1,y2)

# Gamma

alpha1 <- 20
alpha2 <- 10

lam1 <- 5
lam2 <- 10

x1 <- rgamma(n, shape = alpha1, rate = lam1)
x2 <- rgamma(n, shape = alpha2, rate = lam2) 

X11()
plot(x1,x2)

y1 <- x1 + x2
y2 <- x1 - x2

X11()
plot(y1,y2)

# Normal

x1 <- rnorm(n)
x2 <- rnorm(n)

X11()
plot(x1,x2)

y1 <- x1 + x2
y2 <- x1 - x2

X11()
plot(y1,y2)

# Multivariate Normal

mu <- c(50, 100)

VCmatrix <- matrix(c(1,.45,.45,1),ncol=2)

x <- rmvnorm(n, mean = mu, sigma = VCmatrix )

X11()
plot(x[,1],x[,2])

y1 <- x[,1] + x[,2]
y2 <- x[,1] - x[,2]

X11()
plot(y1,y2)

################################################################
# Below is an S-Plus program that plots the joint density of two independent 
# exponential random variables.  What are the values of the rates?  

# Modify the function f to plot the joint density two independent 
# standard normal random variables.  

### Example of a joint density.

x1 <- seq(0,5,0.2)
x2 <- x1
f <- function(x,y){ 2*exp(-x)*exp(-2*y) }
z <- outer(x1,x2,f)

X11()
persp(z)

################################################################
### Example of plotting the joint normal distribution. 

# Compare the perspective plots to the contour plots.

x <- seq(-2,2,0.1)
y <- seq(-2,2,0.1)

X11()
par(mfrow=c(2,2))

bvn <- function(x,y){
    mu1 <- 0
    sigma1 <- 1
    mu2 <- 0
    sigma2 <- 1
    rho <- 0
    ((2*pi*sigma1*sigma2*sqrt(1-rho^2))^(-1))*
            exp( -(2*((1-rho^2))^(-1)) * ( (x-mu1)^2/sigma1^2 + (y-mu2)^2/sigma2^2
            - ( (2*rho*(x-mu1)*(y-mu2) )/ (sigma1*sigma2) ) ) )
}

z1 <- outer(x,y,bvn)
persp(x,y,z1)

bvn <- function(x,y){
    mu1 <- 0
    sigma1 <- 1
    mu2 <- 0
    sigma2 <- 1
    rho <- 0.3
    ((2*pi*sigma1*sigma2*sqrt(1-rho^2))^(-1))*
            exp( -(2*((1-rho^2))^(-1)) * ( (x-mu1)^2/sigma1^2 + (y-mu2)^2/sigma2^2
            - ( (2*rho*(x-mu1)*(y-mu2) )/ (sigma1*sigma2) ) ) )
}

z2 <- outer(x,y,bvn)
persp(x,y,z2)

bvn <- function(x,y){
    mu1 <- 0
    sigma1 <- 1
    mu2 <- 0
    sigma2 <- 1
    rho <- 0.6
    ((2*pi*sigma1*sigma2*sqrt(1-rho^2))^(-1))*
            exp( -(2*((1-rho^2))^(-1)) * ( (x-mu1)^2/sigma1^2 + (y-mu2)^2/sigma2^2
            - ( (2*rho*(x-mu1)*(y-mu2) )/ (sigma1*sigma2) ) ) )
}

z3 <- outer(x,y,bvn)
persp(x,y,z3)

bvn <- function(x,y){
    mu1 <- 0
    sigma1 <- 1
    mu2 <- 0
    sigma2 <- 1
    rho <- -0.9
    ((2*pi*sigma1*sigma2*sqrt(1-rho^2))^(-1))*
            exp( -(2*(1-rho^2)^(-1)) * ( (x-mu1)^2/sigma1^2 + (y-mu2)^2/sigma2^2
            - ( (2*rho*(x-mu1)*(y-mu2) )/ (sigma1*sigma2) ) ) )
}

z4 <- outer(x,y,bvn)
persp(x,y,z4)

X11()
par(mfrow=c(2,2))

contour(x,y,z1)
contour(x,y,z2)
contour(x,y,z3)
contour(x,y,z4)