function s_values=cihm(data,f,g,h,theta,lambda,hpar,m,b)



% GIBBS SAMPLING FOR CONDITIONALLY INDEPENDENT HIERARCHICAL MODEL

% MODEL:

%   	Y(1),...,Y(N) INDEPENDENT f(YI | THETA(I))

%   	THETA(1),...,THETA(N) IID g(THETA | LAMBDA)

%   	LAMBDA DISTRIBUTED h(LAMBDA | HPAR)

% INPUT:

%	FUNCTIONS f,g,h

%	STARTING VALUES theta AND lambda

%   	VECTOR OF HYPERPARAMETERS hpar

%   	NUMBER OF ITERATIONS m TO COLLECT

%	NUMBER OF BURN-IN INTERATIONS b

% OUTPUT:

%   	MATRIX OF SIMULATED VALUES s_values

%   	STORED AS [THETA LAMBDA] WITH m ROWS

% EXAMPLE 

%	(FOR BINOMIAL/BETA EXCHANGEABLE MODEL)

%	s_values=cihm(data,'binomial','beta','hprior',theta,lambda,hpar,1000);

														

K=length(theta); L=length(lambda);

c1=ones(K,1); c2=ones(1,L);	

if nargin==8, b=100; end					%  INITIAL METROPOLIS STEPS

s_values=zeros(b,K+L); 

								%  DEFAULT 100 BURN-IN ITERATIONS

for i=1:(m+b)



	lx=feval(f,data,theta)+feval(g,theta,lambda);			% 

	theta0=theta+randn(K,1).*c1;					%

	ly=feval(f,data,theta0)+feval(g,theta0,lambda);			%  SIMULATES [THETA | LAMBDA]

	pr=exp(ly-lx); in=rand(K,1)<pr;					%

	theta=theta0.*(in==1)+theta.*(in==0);				%



	lx=sum(feval(g,theta,lambda))+feval(h,lambda,hpar);		%

	lambda0=lambda+randn(1,L).*c2;					%

	ly=sum(feval(g,theta,lambda0))+feval(h,lambda0,hpar);		%  SIMULATES [LAMBDA | THETA]

	pr=exp(ly-lx); in=rand<pr;					%

	lambda=lambda0*(in==1)+lambda*(in==0);				%



 	if i<=b        

	 s_values(i,:)=[theta' lambda];					%  COLLECT SIMULATED VARIATES

        else

         s_values(i-b,:)=[theta' lambda];	

        end

	if i/100==fix(i/100),(i-b),end					%  PRINTS CURRENT ITERATION NUMBER

	

	if i==b								%

		s=std(s_values); c1=2*s(1:K)'; c2=s((K+1):(K+2));	%  UPDATES SIZE OF METROPOLIS STEPS

		s_values=zeros(m,K+L); c1,c2				%  AFTER BURN-IN PERIOD

	end								%

end

