
function[deltareturn,accept] = OPdrawdelta(y,x,beta,delta0,d0,D0,tune,Dhat)

% The function samples the cut-point vector using random-walk MH algorithm.
% Output
%--------------------------------------------------------------------------
% deltaReturn : delta from MH algorithm
% accept    : indicator for acceptance of proposed value of delta
% Input
% y         : observed ordinal values, a column vector
% x         : covariate matrix including a column of ones
% beta      : Gibbs draw of beta, a column vector
% delta0    : initial values of cut-points, a column vector
% d0        : prior mean of delta, a column vector
% D0        : prior covariance matrix of delta
% tune      : tuning parameter to adjust MH acceptance rate
% Dhat      : negative inverse Hessian from maximization of log-likelihood
%--------------------------------------------------------------------------

k   = size(delta0,1);

% Proposed value
L      = chol(Dhat,'lower');

delta1 = delta0 + tune*L*normrnd(0,1,k,1);

% Target density

pnum    = -OPnegSumloglike(delta1,y,x,beta) + log(mvnpdf(delta1,d0,D0));

qnum    = -OPnegSumloglike(delta0,y,x,beta) + log(mvnpdf(delta0,d0,D0));

if log(rand) <= pnum-qnum
    deltareturn = delta1;
    accept      = 1;
else
    deltareturn = delta0;
    accept      = 0;
end

end
