
function y = rndtruncNormal(mu,s2,a,b)

% This function draws random numbers from a truncated normal distribution.
% The truncation can be one-sided or two sided and accommodates negative and
% positive infinity. For values which are closer to the mean
% (abs((y-mu)/s)<5), it uses the inverse cdf method. For values which are
% farther away, it uses the accept-reject method using translated
% exponential distribution as the proposal density. See Robert (1995),
% Simulation of Truncated Normal Variables, Statistics and Computing, 5,
% 121-125. See Ivan's Gauss Codes.

% OUTPUT
%--------------------------------------------------------------------------
% y     : draw from a truncated normal distribution
%
% INPUT
%--------------------------------------------------------------------------
% mu    : mean of the untruncated normal distribution
% s2    : variance of the untruncated normal distribution
% a     : left truncation point
% b     : right truncation point
%--------------------------------------------------------------------------

sigma   = sqrt(s2);

if a == -Inf
    if abs((mu-b)/sigma) > 5
        y = 2*mu - rndtruncRob95pos(mu,s2,2*mu-b);
    else
        u     = rand(1,1);
        aStd  = (a-mu)/sigma;    
        bStd  = (b-mu)/sigma;
        wtemp = fast_normcdf(aStd) + u*( fast_normcdf(bStd) - fast_normcdf(aStd));
        y     = mu + sigma.*fast_norminv(wtemp);
    end
%----------
elseif b == Inf
    if abs((mu-b)/sigma) > 5
        y = rndtruncRob95pos(mu,s2,a);
    else
        u     = rand(1,1);
        aStd  = (a-mu)/sigma;    
        bStd  = (b-mu)/sigma;
        wtemp = fast_normcdf(aStd) + u*( fast_normcdf(bStd) - fast_normcdf(aStd));
        y     = mu + sigma.*fast_norminv(wtemp);
    end
%----------
else 
    mua = (mu - a);
    mub = (mu - b);
    % If on the same side of mu and far in the tails
    %-----
    if (mua*mub > 0) && (abs(mua/sigma) > 5) && (abs(mub/sigma) > 5)
    % /Do MODIFIED Robert (1995) -- proposal is truncated exponential/
    % /Original used either uniform proposal or A-R with untruncated
    % proposal/
    if a > mu                % /if density is to the right of mu /
        a1 = (a-mu)/sigma;   % /standardize/
        b1 = (b-mu)/sigma - a1;
    elseif b <= mu           % /if denisty to the left of mu/
        a1 = (mu - b)/sigma;
        b1 = (mu - a)/sigma - a1;
    end
    alstar = (a1 + sqrt(a1^2 + 4))/2; % /optimal inverse scale for exponential proposal/
    % /produce truncation point for uniform u ~ U(b1,1)/ 
    % /leads to sampling from untranslated truncated exponential density (0,b)/
    bL = exp(-alstar*b1);
    
    accept = 0;
    while accept < 1
        u = bL + (1-bL)*rand(1,1);    %/generate from U(bL,1)/
        z = - log(u)/alstar + a1;     % /z ~ rndexpinv(alstar,1)*1{z in (a1,b1)}/
        rho = - ((z-alstar)^2)/2;
        accept = log(rand(1,1)) < rho;
    end
    
    if a > mu                         %/density to the right of mu/
        y = mu + sigma*z;             %/back to TN(mu,s2)/
    elseif b <= mu                    %/density to the left of mu/
        y = mu - sigma*z;
    end
    %-----
    else                 %/if in the stable range of normcdf and norminv/
        u     = rand(1,1);
        aStd  = (a-mu)/sigma;    
        bStd  = (b-mu)/sigma;
        wtemp = fast_normcdf(aStd) + u*( fast_normcdf(bStd) - fast_normcdf(aStd));
        y     = mu + sigma.*fast_norminv(wtemp);
    end
end     


end