%----------------------------------------------------------
% MATLAB for fitting exchangeable model for rat tumour data
%----------------------------------------------------------

% need MATLAB functions betabin, rdisc
%    betabin - function to compute logarithm of posterior density
%    rdisc - function to simulate from a discrete prob. distn.

% sets up grid and computes log posterior on the grid
% graphs posterior by contour plot

 npt=50;
 t1=linspace(-2.3,-1.3,npt);
 t2=linspace(1,5,npt);
 [T1,T2]=meshgrid(t1,t2);
 z=reshape(betabin([T1(:) T2(:)]),npt,npt);
 z=exp(z-max(max(z)));
 figure(1)
 contour(t1,t2,z,.05:.1:.95)
 xlabel('log(alpha/beta)');ylabel('log(alpha+beta)');

% Take a simulated sample of 1000 from this grid

 z=z/sum(sum(z));
 i=rdisc(1000,z(:));
 sim_t1=T1(:); sim_t1=sim_t1(i);
 sim_t2=T2(:); sim_t2=sim_t2(i);
 figure(2)
 plot(sim_t1,sim_t2,'or')
 axis([-2.3 -1.3 1 5])
 xlabel('log(alpha/beta)');ylabel('log(alpha+beta)');

% Transform back to get a simulated sample 
% from posterior of (alpha, beta)

 sim_b=exp(sim_t2)./(1+exp(sim_t1));
 sim_a=exp(sim_t1+sim_t2)./(1+exp(sim_t1));
 figure(3)
 plot(sim_a,sim_b,'ob')
 xlabel('alpha');ylabel('beta')

% Simulate values of (theta(1), ..., theta(71)):
% Plot posterior medians and 95% interval estimates of 
% the theta(i) against the observed proportions y(i)/n(i)

x=[0  0  0  0  0  0  0  0  0  0 ...
 0  0  0  0  1  1  1  1  1  1 ...
 1  1  2  2  2  2  2  2  2  2 ...
 2  1  5  2  5  3  2  7  7  3 ...
 3  5  9 10  4  4  4  4  4  4 ...
 4 10  4  4  4  5 11 12  5  5 ...
 6  5  6  6  6  6 16 15 15  9  4];

n=[ 20 20 20 20 20 20 20 19 19 19 ...
19 18 18 17 20 20 20 20 19 19 ...
18 18 25 24 23 20 20 20 20 20 ...
20 10 49 19 46 27 17 49 47 20 ...
20 16 48 50 20 20 20 20 20 20 ...
20 48 19 19 19 22 46 49 20 20 ...
23 19 22 20 20 20 52 47 46 24 14];

 sim_th=zeros(1000,71);
 for i=1:71, sim_th(:,i)=betarnd(sim_a+x(i),sim_b+n(i)-x(i));i,end
 sim_th=sort(sim_th);
 med=median(sim_th);
 lo=sim_th(25,:);
 hi=sim_th(976,:);
 figure(4)
 errorbar(x./n,med,med-lo,hi-med,'og')
 axis([-.02 .45 -.02 .45]);line([0 .45],[0 .45]);
 xlabel('observed rate y(i)/N(i)')
 ylabel('95% posterior interval for theta(i)')