note10_regression

math 440 linear regression given current values

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Nk (β0 , Σ0 ) Inv-Gamma 2 ν0 σ0 , 22 The joint posterior distribution p(β, σ 2 |y, X ) can be approximated via Gibbs sampling. & % MATH-440 Linear Regression ' $ Given current values {β (t) , σ 2(t) }, new values can be generated by: 1. sampling β (t+1) from its full conditional β |y, X, σ 2(t) ∼ Nk (m(t) , V (t) ), where Slide 31 V (t) = Σ− 1 + 0 1 σ XT X 2(t) −1 m(t) = V (t) Σ−1 β0 + 0 1 σ 2(t) XT y 2. sampling σ 2(t+1) from its full conditional ν σ 2 +S 2(t) + σ 2 |y, X, β (t) ∼ Inv-Gamma ν02 n , 0 0 2 , where S 2(t) = (y − Xβ (t+1) )T (y − Xβ (t+1) ) & % ' $ Let us revisit the bird extinction example. Suppose we take 2 β0 = (0, 0, 0, 0), Σ0 = I , ν0 = 0.2, σ0 = 0.2. Slide 32 b0 = rep(0,k); sigma0 = diag(k); n0=0.2; s0=0.2 T = 10000; nburn = 5000 beta = matrix(NA, T, k); sigma2 = rep(NA, T) beta[1,] = rep(0,k); sigma2[1] = 1 for(t in 2:T) { V = solve(solve(sigma0) + 1/sigma2[(t-1)]*t(x)%*%x) m = V%*%(solve(sigma0)%*%b0 + 1/sigma2[(t-1)]*t(x)%*%y) beta[t,] = rmnorm(1, m, V) S2 = t(y-x%*%beta[t,]) %*% (y-x%*%beta[t,]) sigma2[t] = rigamma(1, (n0+n)/2, (n0*s0 + S2)/2) } & % MATH-440 Linear Regression ' apply(beta, 2, mean) [1] 0.4133903 0.2669139 -0.6282731 Slide 33 $ 0.4977767 apply(beta, 2, quantile, c(0.025, 0.5, 0.975)) [,1] [,2] [,3] [,4] 2.5% 0.007325563 0.1958596 -0.9511456 0.1392730 50% 0.412625341 0.2668820 -0.6270184 0.4985722 97.5% 0.802677165 0.3395624 -0.3071078 0.8464564 quantile(sigma2, c(.025, 0.5, 0.975)) 2.5% 50% 97.5% 0.3015902 0.4289842 0.6267141 & %...
View Full Document

This note was uploaded on 02/22/2013 for the course MATH 440 taught by Professor Tadesse during the Spring '13 term at Georgetown.

Ask a homework question - tutors are online