% the program establishes some data, calls to the Gibbs MCMC
% function for multiple linear regression, and then gets and plots results
% Written by Brice Semmens, 4/10/08
%updated by BXS 7/-2/08

%DATA EXAMPLE -- (taken from course materials by Lawrence Joseph @ McGill)
%Simple linear regression of dependent variable blood pressure (bp) on age and sex.
%Data were simulated from: bp = -25 + 2 * sex + 3 * age, so alpha = -25, beta1 =
%2, beta = 3. Sigma was set to 5. Posterior means should be somewhat near these
%values, but with small sample size, do not expect high accuracy.

Y= [143
132
88
98
177
102
154
83
131
150
131
69
111
114
170
117
96
116
131
158
156
75
111
184
141
182
74
87
183
89];

X=[1	0	59
1	1	52
1	0	37
1	0	40
1	0	67
1	1	43
1	0	61
1	1	34
1	1	51
1	0	58
1	1	54
1	0	31
1	0	49
1	0	45
1	0	66
1	1	48
1	0	41
1	1	47
1	0	53
1	0	62
1	1	60
1	0	33
1	1	44
1	0	70
1	0	56
1	0	69
1	1	35
1	1	36
1	1	68
1	1	38];

%set up priors
priors = zeros(4,2);
priors(1,:) = [.0001,.0001]; %prior for model precision (gam(A,B))
priors(2:4,2) = 10000; %prior variance for the B values
% Now send the data to the Gibbs function:

%set up inits
inits = [100; 1; 1; 1];

% set up model function parmeters
iters = 10000;
burn = 1000;
thin = 10;
tic
posterior=MLR_Gibbs_func(Y,X,inits,priors,burn,thin,iters);
toc
%plot posterior chain
plot(posterior(:,1));