
function [z] = drawlatent(y,x,beta,w,theta,tau2,delta)

% Output
%--------------------------------------------------------------------------
% z         : Gibbs draw of latent response variable, a column vector
% Input
% y         : Ordered discrete values, a column vector
% x         : covariate matrix including a column of ones
% beta      : Gibbs draw of beta, a column vector
% w         : latent weights, a column vector
% theta     : value from normal-exponential mixture of Laplace distribution
% tau       : value from normal-exponential mixture of Laplace distribution
% delta     : cutpoint, modified as delta=log( gam(j) - gam(j-1) )
%--------------------------------------------------------------------------

J          = size(unique(y),1);           % number of outcomes
n          = size(y,1);
z          = zeros(n,1);

% Check function
if J<4
   error('Number of outcomes need to be ge 4 ');
end

expdelta   = exp(delta);                    % for J=4, expdelta is 2x1
gammaCp    = zeros(size(delta,1)+1,1);      % for J=4, gamma =[0 g1 g2]'; gammaCp are the cutpoints.
                                            % Do not use gamma, its a
                                            % function in Matlab
for j = 2:J-1
    gammaCp(j)   = sum(expdelta(1:j-1,1));
end

% Redefining Gamma to include -Inf and Inf
gammaCp       = [-Inf; gammaCp; Inf];

for i =1:n    
    mean    = x(i,:)*beta + theta*w(i);
    sd      = sqrt(tau2*w(i));
   
    temp    = y(i);
    a       = gammaCp(temp);               % left truncation point
    b       = gammaCp(temp+1);             % right truncation point
    %z(i,1)  = rnddtruncnorm(mean,sd,a,b);
    z(i,1) = rndtruncNormal(mean, sd^2,a,b);

end

end
