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



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

% parameters:

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



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

[N,p]=size(y);

mu=mean(y')';					   % STARTING

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

m_al=mean(al); v_al=std(al).^2; v0=1;

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



if nargin<3, ITER=1000; end

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

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

%Mam=zeros(ITER,1); Mav=zeros(ITER,1);



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)+1/v_al);        %

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

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



  v=v_al/N;                                         %

  m=mean(al);                                       % SIMULATE M_AL

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



  a=(v0+N)/2;                                       %

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

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

    

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

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

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

%  Mam(k)=m_al; Mav(k)=v_al;			    %

  Mp1=Mp1+p1; Mp2=Mp2+p2;



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

end

Mp1=Mp1/ITER; Mp2=Mp2/ITER;	

