%
% birats example from WinBUGS
%

dataStruct = struct( ...
    'x', [8.0, 15.0, 22.0, 29.0, 36.0], ...
    'N', 30, ...
    'T', 5, ...
    'Omega', [200, 0; 0, 0.2], ...
    'mean', [0,0], ...
    'prec', [1.0e-6, 0; 0, 1.0e-6], ...
    'Y', [ 151, 199, 246, 283, 320; ...
	   145, 199, 249, 293, 354; ...
	   147, 214, 263, 312, 328; ...
	   155, 200, 237, 272, 297; ...
	   135, 188, 230, 280, 323; ...
	   159, 210, 252, 298, 331; ...
	   141, 189, 231, 275, 305; ...
	   159, 201, 248, 297, 338; ...
	   177, 236, 285, 350, 376; ...
	   134, 182, 220, 260, 296; ...
	   160, 208, 261, 313, 352; ...
	   143, 188, 220, 273, 314; ...
	   154, 200, 244, 289, 325; ...
	   171, 221, 270, 326, 358; ...
	   163, 216, 242, 281, 312; ...
	   160, 207, 248, 288, 324; ...
	   142, 187, 234, 280, 316; ...
	   156, 203, 243, 283, 317; ...
	   157, 212, 259, 307, 336; ...
	   152, 203, 246, 286, 321; ...
	   154, 205, 253, 298, 334; ...
	   139, 190, 225, 267, 302; ...
	   146, 191, 229, 272, 302; ...
	   157, 211, 250, 285, 323; ...
	   132, 185, 237, 286, 331; ...
	   160, 207, 257, 303, 345; ...
	   169, 216, 261, 295, 333; ...
	   157, 205, 248, 289, 316; ...
	   137, 180, 219, 258, 291; ...
	   153, 200, 244, 286, 324 ] );    

init = struct( ...
    'mu_beta', [0,0], ...
    'tauC', 1, ...
    'beta', ...
     [100,6;100,6;100,6;100,6;100,6; ...
	 100,6;100,6;100,6;100,6;100,6; ...
	 100,6;100,6;100,6;100,6;100,6;  ...
	 100,6;100,6;100,6;100,6;100,6; ...
	 100,6;100,6;100,6;100,6;100,6; ...
	 100,6;100,6;100,6;100,6;100,6], ...
    'R', [1,0;0,1]);

clear initStructs
initStructs(1) = init;
initStructs(2) = init;
initStructs(3) = init;

[samples, stats, structArray] = matbugs(dataStruct, ...
		fullfile(pwd, 'birats_model.txt'), ...
		'init', initStructs, ...
		'nChains', 3, ...
		'view', 1, 'nburnin', 1000, 'nsamples', 1000, ...
		'thin', 2, ...
		'monitorParams', { 'beta', 'mu_beta', 'R', 'mu', 'sigma' }, ...
		'Bugdir', 'C:/Program Files/WinBUGS14');

    
% Now for the plots.        


figure; hold on;
colours = 'rgb';
for c=1:3
   plot(samples.mu_beta(c,:,1), colours(c));
end   
title('mu beta(1)');


[N T] = size(dataStruct.Y);
weights = dataStruct.Y(:);
ages    = repmat(dataStruct.x, [N 1]);
ages = ages(:);


figure; hold on;
plot(ages, weights, 'r+');

% Now plot the predicted response for each rat as a function of time.
t = 0:1:40;
for i=1:dataStruct.N
    weight_lin = stats.mean.beta(i, 1) + t * stats.mean.beta(i, 2);
    plot(t, weight_lin, 'g:');
end

% Now plot the predicted response of the population model
weight_lin = stats.mean.mu_beta(1) + t * stats.mean.mu_beta(2);
h=plot(t, weight_lin, 'k'); set(h, 'linewidth', 3)

title('red = data, green = sampled coef, black = mean coef')
xlabel('age')
ylabel('weight')



Y = dataStruct.Y;
x = dataStruct.x;
colors = hsvrand(N);
figure(1);clf
for i=1:N
  %subplot(5,6,i)
  h=plot(x, Y(i,:), 'o-');
  set(h, 'color', colors(i,:))
  hold on
end
xlabel('age')
ylabel('weight')
title('birats: raw data')