assignment

Let us assume that \(\sigma = 1\) is known.

Consider \(\pi(\beta) = \phi(\beta;(0,0), 10I_{2})\)

You already know \(\pi(\beta|y)\)

• Find exact value of \(E(\beta_{0}|y), E(\beta_{1}|y), V(\beta_{0}|y)\) and \(V(\beta_{1}|y)\)

library(MASS)
set.seed(1)
n=100
e <- rnorm(n) # residual
beta0 <- 1; beta1 <- 2; beta <- rbind(beta0,beta1) # parameters beta0 and beta1
x <- cbind(rep(1,n),rnorm(n)) #
y <- x %*% beta + e # linear model equation

v=10 # assume beta~ N((0,0), 10*diag(2)) according to the question, so v = 10
mu_beta <- t(solve(t(x) %*% x + 1/v * diag(2)) %*% (t(x) %*% y))
var_beta <- t(solve(t(x) %*% x + 1/v * diag(2)))
mu_beta;var_beta

##         [,1]     [,2]
## [1,] 1.10766 1.996823

##              [,1]         [,2]
## [1,] 0.0100056979 0.0004153511
## [2,] 0.0004153511 0.0109967641

• Use Monte Carlo Method (size = 100, 1000, 5000)

E_beta_mc <- array(NA,dim = c(3,2))
V_beta_mc <- E_beta_mc
size <- c(100, 1000, 5000)
for (i in 1:3) {
# sample from \pi(\beta|y)
beta_mc <- as.data.frame(mvrnorm(n=size[i], mu_beta, var_beta)) 
names(beta_mc) <-c('beta0','beta1') 
E_beta_mc[i,] <- colMeans(beta_mc)
V_beta_mc[i,] <- apply(beta_mc,2,var)
}
E_beta_mc

##          [,1]     [,2]
## [1,] 1.103916 2.001518
## [2,] 1.107062 1.992977
## [3,] 1.108487 1.995926

V_beta_mc

##             [,1]       [,2]
## [1,] 0.009187138 0.01242899
## [2,] 0.010760386 0.01235912
## [3,] 0.009881541 0.01108607

• Importance Sampling

Use \(\phi(\beta_{0} |1, 0.5),\phi(\beta_{1} |2, 0.5)\)

E_h_beta <- array(NA, dim = c(3,2))
V_h_beta=E_h_beta
for (i in 1:3) {
# g_beta0 represents beta0 sampled from N(1,0.5)
g_beta0 <- rnorm(size[i], mean = 1, sd = sqrt(0.5)) 
# g_beta1 represents beta1 sampled from N(2,0.5)
g_beta1 <- rnorm(size[i], mean = 2, sd = sqrt(0.5))
# density function value from g(\beta0)
density_beta0 <- dnorm(g_beta0, mean = 1, sd = sqrt(0.5)) 
# density function value from g(\beta1)
density_beta1 <- dnorm(g_beta1, mean = 2, sd = sqrt(0.5)) 
# density function value from prior dist'n
pi_beta0 <- dnorm(g_beta0, mean = mu_beta[1], sd = sqrt(var_beta[1,1])) 
# density function value from prior dist'n
pi_beta1 <- dnorm(g_beta1, mean = mu_beta[2], sd = sqrt(var_beta[2,2])) 

h_beta0 <- g_beta0*pi_beta0/density_beta0
h_beta1 <- g_beta1*pi_beta1/density_beta1

E_h_beta0 <- mean(h_beta0)
V_h_beta0 <- var(h_beta0)
E_h_beta1 <- mean(h_beta1)
V_h_beta1 <- var(h_beta0)
E_h_beta[i,] <- c(E_h_beta0,E_h_beta1)
V_h_beta[i,] <- c(V_h_beta0,V_h_beta1)
}
E_h_beta

##          [,1]     [,2]
## [1,] 1.050015 2.062285
## [2,] 1.132268 1.952093
## [3,] 1.126902 2.047798

V_h_beta

##          [,1]     [,2]
## [1,] 4.415393 4.415393
## [2,] 5.219405 5.219405
## [3,] 5.097783 5.097783