function [Zij,Z,S,Cats,B,Sr,accept]=sampReg(N,X,bZ,sampSize,alpha,lambda)
% sampleReg returns a sample from the posterior on the
%   rater-item values (Zij), the latent item traits, Z, the rater
%   variances S, and the category-cutoffs Cats.   It also generates
%   a sample from the regression vector of Z on X, the design matrix
%   IT IS ASSUMED THAT THE MULTINOMIAL DENOMINATORS ARE ALL 1.
%   
%   N:		the data matrix, with "items" rows and "raters" columns.
%               Entries are labeled 1 to C, and it is assumed that each
%		rater rated each item.     
%   X:          Design matrix, (items x nparm)
%   bZ:	        estimate of mean latent trait vector, (items x 1)
%   sampSize:	desired number of MCMC iterates. 
%   alpha,lambda: prior parameters on rater variances 
%   Zij:	sampSize x items x raters array containing sample
%               rater-item values z_ij.
%   Z:          sampSize x items array containing estimates z_i's
%   S           sampSize x rater array wit sampled variances
%   Cats        sampSize x C+1 x rater array with sampled category
%               cutoffs.  For simplicity (,1,) element is -10,
%               (,C+1,1)=10.  (i,j,r) is the ith sampled value of the
%               lower cutoff for the j'th category (upper cutoff for
%               the (j-1) category) for the r'th rater
%   B:	        Sampled values of the regression parameter
%   Sr:	        Sampled values of the regression variance
%   
%   accept:	acceptance ratio for category cutoffs

%   Matlab initialize vectors

    [items raters] = size(N);
    [items nparm] = size(X);
    C = max(max(N));
    Zij(1:sampSize,1:items,1:raters) = 0;
    Z = zeros(sampSize,items);
    S = zeros(sampSize,raters);
    Cats(1:sampSize,1:(C+1),1:raters) = 0;
    B = zeros(sampSize,nparm);
    Sr = zeros(sampSize,1);

    H = X*inv(X'*X)*X';
    Hperp = eye(size(H))-H;
    H0 = inv(X'*X)*X';
    XtX = X'*X;
    iXtX = inv(XtX);
    B(1,:) = (H0*bZ)';

    Sr(1) = 1;

%   "Sample" first value of vectors

    muN = mean(N);
    stdN = sqrt(var(N));
    Zij(1,:,:) = (N-ones(items,1)*muN)./(ones(items,1)*stdN);
    temp = squeeze(Zij(1,:,:));
    Z(1,:) = mean(temp');

    Cats(:,1,:) = -10;
    Cats(:,(C+1),:) = 10;
    for i=2:C
	Ind = ones(size(N));
        j = find(N>=i);
	Ind(j) = 0;
        Cats(1,i,:) = Phiinv((sum(Ind)+.5)/(items+.501));
    end

    tZij = squeeze(Zij(1,:,:));
    tZ   = squeeze(Z(1,:))'; 
    ss   = sum( (tZij-tZ*ones(1,raters)).^2 )';
    temp = (rgamma( 0.5*items + alpha*ones(raters,1), 0.5*ss +lambda)).^(-1);
    S(1,:) = temp';

    sigma = sqrt(S);
    
    sm = 0.4/C;
    g = zeros(C+1,1);
    oldg = zeros(C+1,1);
    accept = 0;


%   Begin real sampling loop;

    for i=2:sampSize

%     1. Sample latent vector Zij

	  for j=1:raters
	     upper = Cats((i-1),(N(:,j)+1),j)';
	     lower = Cats((i-1),N(:,j),j)';
             Zij(i,:,j) = truncNorm(Z((i-1),:)',sigma((i-1),j),lower,upper); 
	  end

%     2. Sample gamma
          
	  for j=1:raters      
%	      a) Put proposal gamma in g; the old gamma in oldg
		oldg = Cats((i-1),:,j)'; 
                g = oldg;
                for k=2:C 
		   g(k) = truncNorm(oldg(k),sm,g(k-1),oldg(k+1));
	        end
	    

%	      b) Calculate acceptance ratio R
		%   adjust R for proposal density truncation
                R = 1;
		for k=2:C
		   R = R * ( Phi((oldg(k+1)-oldg(k))/sm) -  ...
                         Phi((g(k-1)-oldg(k))/sm) ) / ... 
                       ( Phi((g(k+1)-g(k))/sm) - ...
                         Phi((oldg(k-1)-g(k))/sm) );
	        end
	        %  multiply in likelihood

                phiYnew = Phi( (g(N(:,j)+1) - Z((i-1),:)')/sigma((i-1),j) );
                phiYold = Phi( (oldg(N(:,j)+1) - Z((i-1),:)')/sigma((i-1),j) );
		phiYm1new=Phi( (g(N(:,j)) - Z((i-1),:)')/sigma((i-1),j));
	        phiYm1old=Phi( (oldg(N(:,j)) -Z((i-1),:)')/sigma((i-1),j));
	        R = R*prod( (phiYnew-phiYm1new)./(phiYold-phiYm1old));

%             c) Accept/reject
                % accept/reject
                if rand < R
                  Cats(i,:,j) = g';
                  accept = accept+1;
		else
		  Cats(i,:,j) = oldg';
                end

	    end

  
%     3. Sample Z given Zij
                b = sum((squeeze(Zij(i,:,:))./(ones(items,1)*S((i-1),:)))' )';
                a = 1+sum((ones(items,raters)./(ones(items,1)*S((i-1),:)))' )';
%		Z(i,:) = (b ./ a)'+randn(1,items)./sqrt(a');
%
%          Do MH accept step to account for regression parameter
%
		temp = (b ./ a)'+randn(1,items)./sqrt(a'); %new value
                t0 = (temp'-X*B((i-1),:)')'*(temp'-X*B((i-1),:)');
                told = (Z(i,:)'-X*B((i-1),:)')'*(Z(i,:)'-X*B((i-1),:)');
	        rssNew = temp*Hperp*temp';
                rssOld = Z(i,:)*Hperp*Z(i,:)';
		R = exp( -( t0 - told )/(2*Sr(i-1)) + ...
                         (0.5*(items-nparm)+alpha)*log((rssNew+2*lambda) ...
                                  /  (rssOld+2.0*lambda)) );
	        
	        if R>rand
                   Z(i,:)  = temp;
                else
                   Z(i,:) = Z((i-1),:);
                   rssNew = rssOld;
                end

%     4. Sample S
	        tZij = squeeze(Zij(i,:,:));
		tZ   = squeeze(Z(i,:))'; 
		ss   = sum( (tZij-tZ*ones(1,raters)).^2 )';
		temp = (rgamma( 0.5*items + alpha*ones(raters,1), ...
			0.5*ss +lambda)).^(-1);
		S(i,:) = temp';
		sigma = sqrt(S);

%     5. Sample B
                B(i,:) = rMultiNorm(H0*Z(i,:)',Sr(i-1)*iXtX)';

%     6. Sample Sr
                temp = (rgamma( 0.5*(items-nparm)+alpha, ...
                         0.5*rssNew+lambda) )^(-1);
                Sr(i) = temp;

  end

  accept = accept/((sampSize-1)*raters);

    function val=Phi(x)
val=.5*(1+erf(x/sqrt(2)));

function val=Phiinv(x)
val=sqrt(2)*erfinv(2*x-1);

function gam = rgamma(alpha,lambda)
% Generates random gamma deviates from density
%       lambda^alpha g^alpha-1 exp(-lambda*x)/Gamma(alpha)
%       E(gam) = alpha/lambda   Var(gam)=a/lambda^2
%
%  Begin by generating g~gamma(alpha,1)
%
  e = exp(1);
  gam = -ones(size(alpha));
  if alpha<1
     for i=1:length(alpha)
        g = -1; 
        a = alpha(i);
        aa = (a+e)/e;
        while g<0
          r1 = rand;
          r2 = rand;
          if r1>1/aa
             w = -log(aa*(1-r1)/a);
             if r2 < w^(a-1)
                g = w;
             end
          else
             w = (aa*r1)^(1/a);
             if r2<exp(-w)
                g = w;
             end
          end
        end
        gam(i) = g;
     end
  elseif alpha>1
     for i=1:length(alpha)
        a = alpha(i)-1;
        b = (alpha(i)-1/(6*alpha(i)))/a; 
        m = 2/a;
        d = m+2;
        g = -1;
        while g<0
           x = rand;
           y = rand;
           v = b*y/x;
           if (m*x-d+v+1/v) <= 0
                g = a*v;
           elseif (m*log(x)-log(v)+v-1) <= 0
                g = a*v;
           end
        end
        gam(i) = g;
      end
  else
        gam = -log(rand(size(alpha)));
  end
  gam = gam ./ lambda;
         
function val = truncNorm(mu,std,lower,upper)
% truncNorm returns a sample vector of normal deviates with means mu
%     standard deviation std, truncated to the intervals
%     (lower,upper).
%

% Calculate bounds on probabilities
  lowerProb = Phi((lower-mu)./std);
  upperProb = Phi((upper-mu)./std);

% Draw uniform from within (lowerProb,upperProb)
  u = lowerProb+(upperProb-lowerProb).*rand(size(mu));

% Find needed quantiles
  val = mu + Phiinv(u).*std;

  function rnorm = rMultiNorm(mu,cov)
% rMultiNorm generates a random multivariate normal vector.
%
randn(size(mu));
rnorm = mu +chol(cov)'*randn(size(mu));