forked from burakbayramli/books
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ch11q3.m
52 lines (45 loc) · 1.38 KB
/
ch11q3.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
%Matlab code for Question 3 of Chapter 11
%Load in the data you generated in part a)-- relevant code to be found in ch11q3a.m
load ch2q2.out;
n=size(ch2q2,1);
y=ch2q2(:,1);
x=ch2q2(:,2:3);
%x=ch2q2(:,2);
k=size(x,2);
%Hyperparameters for natural conjugate prior
%I am using the suffix "0" to denote prior hyperparameters
v0=1;
s2inv0=1;
s02=1/s2inv0;
b0=zeros(k,1);
b0(k,1)=1;
c=1;
capv0=c*eye(k);
%Parameters for Normal-Gamma posterior
%First let's do OLS since some posterior results written in terms of OLS
%Ordinary least squares quantities
xsquare=x'*x;
bols = inv(xsquare)*x'*y;
s2 = (y-x*bols)'*(y-x*bols)/(n-k);
v=n-k;
%Now use OLS in the formulae for the posterior
v1=v0+n;
capv0inv = inv(capv0);
capv1inv = capv0inv+ xsquare;
capv1=inv(capv1inv);
b1 = capv1*(capv0inv*b0 + xsquare*bols);
v1s12 = v0*s02 + v*s2 + (bols-b0)'*inv(capv0 + inv(xsquare))*(bols-b0);
s12 = v1s12/v1;
bcov = capv1*v1s12/(v1-k);
bsd=zeros(k,1);
for i = 1:k
bsd(i,1)=sqrt(bcov(i,i));
end
%Now we can print out the posterior mean and standard deviation of beta
'posterior mean and standard deviation of beta'
[b1 bsd]
%log of marginal likelihood for the model if prior is informative
intcon=gammaln(.5*v1) + .5*v0*log(v0*s02)- gammaln(.5*v0) -.5*n*log(pi);
lmarglik=intcon + .5*log(det(capv1)/det(capv0)) - .5*v1*log(v1s12);
'log of marginal likelihood'
lmarglik