# R Code 10, Blocked Sampling

Consider the linear model
Y ~ X theta + epsilon

Y: column vector of N observations
X: design matrix, N by R
theta: R parameters
epsilon: column vector of zero-center iid normals with variance sigma^2

It is sometimes very convenient to update all elements of theta at once with a multivariate draw.

Prior on theta: theta ~ N(mu_0,Sigma0)
mu0: column vector of R means
Sigma0: R-by-R covariance matrix.

f(theta|Y) ~ L(Y|theta) X f(theta)

The derivation requires factoring a matrix and is a bit cumbersome. Here are the very important results.

theta|Y ~ Normal(vc,v),

where v is a R-by-R covariance matrix
v=(X'X/sigma^2+B)^{-1}
and B=Sigma0^{-1}

and where c is the column vector, c=X'y/sigma^2+Bmu_0

################################################

##########################################

library(MASS)
library(MCMCpack)

#Set True Values
set.seed(123)
I=50
J=5
N=I*J
vTest=rnorm(I)
cov=(vTest-mean(vTest))/sd(vTest)
t.beta=.2
t.mu=1
t.alpha=rnorm(I,0,.2)
t.theta=c(t.mu,t.alpha,t.beta)
t.sig2=(.2)^2

###########################
#Design Matrix, fullest version
sub=rep(1:I,each=J)
X=matrix(0,nrow=N,ncol=I+2)
for (r in 1:N)
{
X[r,1]=1
X[r,1+sub[r]]=1
X[r,I+2]=cov[sub[r]]
}
image(X,col=grey((0:256)/256))

rt=rnorm(N,X%*%t.theta,sqrt(t.sig2))
obsSubMean=tapply(rt,sub,mean)
plot(cov,obsSubMean)

###########Bayesian Analysis 1, Regression
# y_i = mu + beta*cov+noise
M=2000

Xr=X[,c(1,I+2)]
theta=matrix(nrow=M,ncol=2)
sig2=1:M
theta[1,]=c(.5,0) #start values
sig2[1]=.01

B=diag(c(.1,.1))
mu0=c(0,0)
a=1
b=.1

XtX.r=t(Xr)%*%Xr
Xty.r=t(Xr)%*%rt
Bmu=B%*%mu0

for (m in 2:M)
{
# sample (theta|sig2,y)
v=solve(XtX.r/sig2[m-1]+B)
c=Xty.r/sig2[m-1]+Bmu
theta[m,]=mvrnorm(1,v%*%c,v)
# sample (sig2|theta,y)
s2.scale=sum((rt-Xr%*%theta[m,])^2)/2+b
sig2[m]=rinvgamma(1,shape=a+N/2,scale=s2.scale)
}

############Bayesian Analysis 2, Shrinkage to Grand Mean, Not Preferred
# y_i = mu+ alpha_{j_i}+noise

Xa=X[,1:(I+1)]
theta=matrix(nrow=M,ncol=(I+1))
sig2=1:M
sig2G=1:M
theta[1,]=c(.5,rep(0,I)) #start values
sig2[1]=.01
sig2G[1]=.01
mu0=rep(0,(I+1))
a=1
b=.1
XtX.a=t(Xa)%*%Xa
Xty.a=t(Xa)%*%rt
for (m in 2:M)
{
B=diag(c(.1,rep(1/sig2G[m-1],I)))
v=solve(XtX.a/sig2[m-1]+B)
c=Xty.a/sig2[m-1]
theta[m,]=mvrnorm(1,v%*%c,v)
s2.scale=sum((rt-Xa%*%theta[m,])^2)/2+b
sig2[m]=rinvgamma(1,shape=a+N/2,scale=s2.scale)
s2G.scale=sum(theta[m,2:(I+1)]^2)/2+b
sig2G[m]=rinvgamma(1,shape=a+I/2,scale=s2G.scale)
}

plot(theta[,1])
post.mean=apply(theta[200:M,],2,mean)
post.sub=post.mean[2:(I+1)]
plot(t.theta[2:(I+1)],post.sub)
plot(obsSubMean-mean(obsSubMean),post.sub)
abline(0,1)
plot(cov,post.sub)
mean(sig2G[200:M]) #unaccounted subject variance

############################################################## #Analysis #2 Build in the Regression Effect into model too
# y_i = mu+alpha_{j_i} + beta*cov_i+epsilon_i
theta=matrix(nrow=M,ncol=(I+2))
sig2=1:M
sig2G=1:M
theta[1,]=c(.5,rep(0,I),0) #start values
sig2[1]=.01
sig2G[1]=.01
a=1
b=.1
XtX=t(X)%*%X
Xty=t(X)%*%rt
for (m in 2:M)
{
B=diag(c(.1,rep(1/sig2G[m-1],I),1))
v=solve(XtX/sig2[m-1]+B)
c=Xty/sig2[m-1]
theta[m,]=mvrnorm(1,v%*%c,v)
s2.scale=sum((rt-X%*%theta[m,])^2)/2+b
sig2[m]=rinvgamma(1,shape=a+N/2,scale=s2.scale)
s2G.scale=sum(theta[m,2:(I+1)]^2)/2+b
sig2G[m]=rinvgamma(1,shape=a+I/2,scale=s2G.scale)
}

plot(theta[,1])
plot(theta[,(I+2)])
post.mean=apply(theta[200:M,],2,mean)
post.sub=post.mean[2:(I+1)]
plot(t.theta[2:(I+1)],post.sub)
abline(0,1)
plot(obsSubMean-mean(obsSubMean),post.sub)
abline(0,1)
plot(cov,post.sub)
mean(sig2G[200:M]) #unaccounted subject variance