% MATH 648 - hierarchical modeling hw - exercise 2
%
% fitting a binomial/beta exchangeable model
%    1. y_i are independent binomial(n_i, p_i)
%    2. p_i are iid beta(Km, K(1-m))
%    3. (K, m) have joint prior 1/(m(1-m)) 1/(K+1)^2
%
% first setup the data in Matlab
% y is number of hits for all players, n are the corresponding number
% of bats

y=[12 7 14 9 18 16 15 13 8 10 17 10 11 10 14 11 10 10]';
n=45*ones(18,1);
data=[y n];

% we use the function betabin.m which contains the definition of the
% log posterior of the hyperparameters
% t1 = log(m/(1-m)), t2 = log(K)

% first use function laplace to find posterior mode and variance-cov matrix

[mu,v]=laplace('betabin',[-1 2],10,data);

% display mode and associated standard deviations

mu, sqrt(diag(v))

ITER=1000; % number of simulation iterations
SIMDATA=zeros(ITER,20); % stores simulated values of m,K,p1,...,p18

for k=1:ITER
	
	t1t2=rbivnorm(mu,v);         % simulates values of logit(m), log(K)
	t1=t1t2(1); t2=t1t2(2);      % from approximate bivariate normal dist.
	m=exp(t1)/(1+exp(t1)); K=exp(t2);

	p=zeros(1,18);          
	for j=1:18                   % simulates pj from indep. beta distn's
		p(j)=betarnd(K*m+y(j),K*(1-m)+45-y(j));
	end

	SIMDATA(k,:)=[m K p]	
	
end