Exercise¶

1. Use $\beta_0=1$, $\sigma^2_{\beta}=0.1$ and $\sigma^2_e = 1.0$ to generate a data set with 10 observations from model (2) with $k=15$ covariates.

2. Setup and solve the mixed model equations given by (3).

3. Sample the elements of $\boldsymbol{\beta}$ using Gibbs.

4. Compute the posterior mean of $\boldsymbol{\beta}$ from the samples and compare with the mixed model solutions.

5. Compute the posterior covariance matrix from the sampled values. Compare results with inverse of the mixed-model coefficient matrix.

Simulation of Data¶

In [12]:
using Distributions
using StatsBase

In [13]:
n = 10 #number of observations
k = 15    #number of covariates

varBeta = 0.1
varRes  = 1.0

x = sample([0,1,2],(n,k))
X = [ones(n) x]

betaTrue=[1, randn(k)*sqrt(varBeta)]
y = X*betaTrue+ randn(n)*sqrt(varRes);


MME¶

In [14]:
λ=varRes/varBeta
d=eye(k+1)
d[1,1]=0
lhs=X'X+d*λ
rhs=X'y
ilhs = inv(lhs)
sol=ilhs*rhs

Out[14]:
16-element Array{Float64,1}:
0.739611
-0.0985318
-0.0733791
0.244959
0.171841
0.0126194
0.0451185
-0.0675262
-0.0307764
-0.181414
-0.193714
0.255449
-0.194699
-0.183016
-0.113338
0.153745 
In [15]:
ilhs*varRes

Out[15]:
16x16 Array{Float64,2}:
1.477      -0.0714228    -0.117995     …  -0.101212     -0.05469
-0.0714228   0.0758958     0.00701886       0.00749083    0.0067372
-0.117995    0.00701886    0.0702412        0.000819962  -0.00193193
-0.0733976   0.00211372    0.00179398       0.00259769   -0.00573441
-0.0891733   0.00406518    0.00509045       0.000641145  -0.00204173
-0.127917    0.000686534   0.0153226    …   0.00329912   -0.0067671
-0.0275212  -0.000976669   0.0147982       -0.0155726     0.00236965
-0.107998   -0.00026764    0.0087143        0.006673      0.00504872
-0.0340253  -0.00482864    0.00436116       0.00322588    0.00317027
-0.0642605  -0.0210254     0.0034216       -0.000252368  -0.00144475
-0.0954171   0.00105101    0.00223758   …  -0.00765011    0.0131453
-0.126453    0.00377223   -0.00606176       0.0139139    -0.00279948
-0.0945257  -0.00416139   -0.00977513      -0.00482177    0.00724818
-0.0921117   0.000430086  -0.00751776      -0.00323865    0.00813793
-0.101212    0.00749083    0.000819962      0.073409     -0.00654971
-0.05469     0.0067372    -0.00193193   …  -0.00654971    0.0779866 

MCMC:¶

In [19]:
# loop for Gibbs sampler
niter  = 100000 # number of samples
burnin = 1000
b     = zeros(1+k)
meanB = zeros(1+k)
vB    = zeros(1+k,1+k)
meansigma2 = 0
nu = 4
sigma2 = 2
covariates = length(b) - 1 #Number of covariates
b_mcmc = zeros(niter, length(b))

for iter = 1:niter

# sampling intercept
w = y - X*b + X[:,1]*b[1]
x = X[:,1]
xpxi = 1/(x'x)[1]
bHat = (xpxi*x'w)[1]
b[1] = rand(Normal(bHat, sqrt(xpxi))) # using residual var = 1

#sample covariates
for j = 2:length(b)
w = y - X*b + X[:,j]*b[j]
x = X[:,j]
xpxi = 1/((x'x)[1]+λ)
bHat = (xpxi*x'w)[1]
b[j] = rand(Normal(bHat, sqrt(xpxi)))
end
if iter >= burnin
b_mcmc[iter,:] = b
meanB += b
vB    += b*b'
end
end


Posterior Mean of $\beta$ Estimated from MCMC Samples¶

In [20]:
meanB /= (niter-burnin)

Out[20]:
16-element Array{Float64,1}:
1.20281
-0.0941883
-0.143509
0.314295
0.222457
-0.108558
0.0911845
-0.149743
-0.00393558
-0.263458
-0.296184
0.348927
-0.210706
-0.341597
-0.152182
0.179394  

Posterior Covariance of $\beta$ Estimated from MCMC Samples¶

In [9]:
vB/(niter-burnin) - meanB*meanB'

Out[9]:
16x16 Array{Float64,2}:
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0

Unknown variance components¶

In [29]:
k=15
niter  = 100000 # number of samples
burnin = 1000
b     = zeros(1+k)
meanB = zeros(1+k)
vB    = zeros(1+k,1+k)
meansigma2 = 0
nu = 4
sigma2 = 2
covariates = length(b) - 1 #Number of covariates

dfEffectVar=4
scaleVar=0.1
chi2=Chisq(dfEffectVar+covariates)

sampleVar = zeros(niter)

for iter = 1:niter

# sampling intercept
w = y - X*b + X[:,1]*b[1]
x = X[:,1]
xpxi = 1/(x'x)[1]
bHat = (xpxi*x'w)[1]
b[1] = rand(Normal(bHat, sqrt(xpxi))) # using residual var = 1

#sample covariates
λ=varRes/varBeta
for j = 2:length(b)
w = y - X*b + X[:,j]*b[j]
x = X[:,j]
xpxi = 1/((x'x)[1]+λ)
bHat = (xpxi*x'w)[1]
b[j] = rand(Normal(bHat, sqrt(xpxi)))
end

varBeta = (scaleVar*dfEffectVar + dot(b[2:end],b[2:end]))/rand(chi2)

if iter >= burnin
meanB += b
vB    += b*b'
end

sampleVar[iter]=varBeta
end

In [34]:
using Gadfly
plot(x=sampleVar,Geom.histogram,Guide.title("Posterior distribution of varBeta"),
Guide.ylabel("Frequency"),
Guide.xlabel("varBeta"))

Out[34]: