
function [w] = drawW(z,x, beta, tau2, theta, lambda)

% Output
%--------------------------------------------------------------------------
% w        : latent weights drawn from GIG distribution, a column vector
% Input
% z         : latent response variable, a column vector
% x         : covariate matrix including a column of ones
% beta      : Gibbs draw of beta, a column vector
% tau2      : parameter from normal-exponential mixture of Laplace dist.
% theta     : parameter from normal-exponential mixture of Laplace dist.
% lambda    : Index parameter of GIG distribution
%--------------------------------------------------------------------------
[n,~]       = size(x);

tildeEta2 = (theta^2)/(tau2) + 2;    % tildeEta2 = tildegamma2, notation in J=3
tildeLambda2 = zeros(n,1);                       % tildeLambda2 = tildedelta2, notation in J=3  

w            = zeros(n,1);
for j = 1:n
    tildeLambda2(j,1)= ( (z(j) - x(j,:)*beta)^2 )/(tau2);
    % Patch
    if  tildeLambda2(j,1) < 10^-8
        tildeLambda2(j,1) = 10^-8;
    end
    % w(j,1)  = randraw('gig',[lambda, tildeLambda2(j,1), tildeEta2],1);
    w(j,1) = gigrnd(lambda, tildeEta2, tildeLambda2(j,1), 1);
    
end
%--------------------------------------------------------------------------  

end