function [Mp1,Mp2]=Tgibbs(y,s2,ITER)



% DATA:  y:  N x 2, sampling variances s2: N x 2

% parameters:

%     mu (N x 1), al (N x 1), m_al, v_al, lam (N x 1)

% hyperparameters:  df, v0



% SETUP ----------------------------------

[N,p]=size(y); df=4; a0=53.2; b0=0.810;

% starting estimates

mu=mean(y')';

res=y-mu*[1 1]; al=res(:,1); 

m_al=mean(al); v_al=std(al).^2; lam=ones(N,1); 

% ----------------------------------------



if nargin<3, ITER=1000; end

Mp1=zeros(ITER,N); Mp2=zeros(ITER,N); 



for k=1:ITER

  z=y+al*[-1 1];														                      %

  v=(1)./((1)./s2(:,1)+(1)./s2(:,2));               %

  m=v.*(z(:,1)./s2(:,1)+z(:,2)./s2(:,2));           % SIMULATE MU                                   %

  mu=m+randn(N,1).*sqrt(v);                         %



  w=y-mu*[1 1];														                           %

  v=(1)./((1)./s2(:,1)+(1)./s2(:,2)+lam/v_al);          %

  m=v.*(w(:,1)./s2(:,1)-w(:,2)./s2(:,2)+lam*m_al/v_al); % SIMULATE ALPHA

  al=m+randn(N,1).*sqrt(v);                             %



  v=v_al/sum(lam);                                  %

  m=sum(lam.*al)/sum(lam);                          % SIMULATE M_AL

  m_al=m+randn*sqrt(v);                             %



  a=N/2+a0;                                         %

  b=sum(lam.*(al-m_al).^2)/2+b0;                    % SIMULATE V_AL

  v_al=1/gamrnd(a,1/b);                             %



  a=(df+1)/2;                                       %

  b=(df+(al-m_al).^2/v_al)/2;                       % SIMULATE LAM

  lam=rgam(N,a)./b;                                 % 



  p1=exp(mu+al)./(1+exp(mu+al));

  p2=exp(mu-al)./(1+exp(mu-al));



  Mp1(k,:)=p1'; Mp2(k,:)=p2';



  if k/50==fix(k/50),k,end

end





	

