% 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));