# R code for HW 4, spring 2014

library(mvtnorm)

# Problem 4:

# y values for data set of California oranges:

y.CA <- c(12.6,10.2,9.2,16.4,12.6,12.6,14.3,14.0,10.3,7.7,13.0,9.8,23.6,29.9,20.0,
24.5,30.3,24.1,21.2,18.2,40.7,44.3,41.4,35.3)

# x values for data set of California oranges:

x.CA <- c(0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,2,2,2,2)


# y values for data set of Florida oranges:

y.FL <- c(3.7,3.8,7.1,6.5,8.5,4.0,10.8,5.2,9.7,0.8,10.2,3.0,6.5,20.6,17.2,
12.8,18.5,15.6,17.3,16.9,16.3,18.8,11.0,31.8,28.6,34.4,29.8,24.3,27.4,25.5)

# x values for data set of Florida oranges:

x.FL <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2)

# California:

y <- y.CA

X <- as.matrix( cbind(rep(1, length(x.CA)), x.CA) )

### Noninformative Prior Analysis:

bhat.CA <- solve(t(X) %*% X) %*% t(X) %*% y
sig2hat.CA <- t(y - X %*% bhat.CA) %*% (y - X %*% bhat.CA) / (nrow(X) - ncol(X))
my.cov.mat.CA <- solve(t(X) %*% X) * ((nrow(X)-ncol(X))*sig2hat.CA / (nrow(X) - ncol(X) - 2))[1,1]

CA.draws <- matrix(bhat.CA,nr=10000,nc=2,byrow=T) + rmvt(n=10000, sigma = my.cov.mat.CA, df = nrow(X) - ncol(X) )

#################################################################
# Florida:

y <- y.FL

X <- as.matrix( cbind(rep(1, length(x.FL)), x.FL) )


### Noninformative Prior Analysis:

bhat.FL <- solve(t(X) %*% X) %*% t(X) %*% y
sig2hat.FL <- t(y - X %*% bhat.FL) %*% (y - X %*% bhat.FL) / (nrow(X) - ncol(X))
my.cov.mat.FL <- solve(t(X) %*% X) * ((nrow(X)-ncol(X))*sig2hat.FL / (nrow(X) - ncol(X) - 2))[1,1]


FL.draws <- matrix(bhat.FL,nr=10000,nc=2,byrow=T) + rmvt(n=10000, sigma = my.cov.mat.FL, df = nrow(X) - ncol(X) )

# You will have to write some code to do part (c) that will use one or more
# of the following randomly generated quantities:

Calif.intercepts <- CA.draws[,1]
Calif.slopes <- CA.draws[,2]
Florida.intercepts <- FL.draws[,1]
Florida.slopes <- FL.draws[,2]

# You will have to write some code to do Part (e) that will also use the randomly generated 
# quantities above.