function y=exprej(a,b)
%Exprej    Sample from the truncated normal using the exponential rejection
%          sampling.



if b==inf
    v=rand;
    u=rand;
    x=-1/a*log(v)+a;
    
    while (u>exp(-1/2*(x-a)^2))
          v=rand;
          x=-1/a*log(v)+a;
          u=rand;
    end

else
    u=rand;
    x=texp(a,b);
    
while (u>exp(-1/2*(x-a)^2)) 
    
    u=rand; 
    x=texp(a,b);   
    
    end

end
    y=x;