new;

library pgraph;graphset;

format /rd 10,8;

rndseed 1166;

timedt;
n = 5000;
x = zeros(n,1);

/*Generate random variates from unknown marginal*/
/*Based on f(x|y) = ye^(-yx)                    */

cand0  = 0;
a = 2;

/*Generate random variates from a Laplace distribution*/
/*Use Normal Distribution to Augment MH Algorithm     */

for j (1, n, 1);

    /*Step 1--Generate a second candidate from f(y)*/

    jj = 0;

    do until jj eq 1;

        u_i = rndu(1,1);
        z_i = rndn(1,1);

        stat = .5*exp(-abs(z_i))*inv(a*pdfn(z_i));

        if stat .ge u_i;

            cand = z_i;
            jj = 1;

        endif;

    endo;

    Cx = (.5*exp(-abs(cand0)) .le a*pdfn(cand0));
    Cy = (.5*exp(-abs(cand)) .le a*pdfn(cand));

    if Cx;

        x[j]  = cand;
        cand0 = cand;

    elseif Cy;

        u_i = rndu(1,1);

        if u_i .le (a*pdfn(cand0)*inv(.5*exp(-abs(cand0))));

            x[j]  = cand;
            cand0 = cand;

        else;

            x[j]  = cand0;
            cand0 = cand0;

        endif;

    else;

        u_i = rndu(1,1);

        stat = .5*exp(-abs(cand))*pdfn(cand0)/(.5*exp(-abs(cand0))*pdfn(cand));
        
        if stat .le 1;

            stat = stat;

        else;

            stat = 1;

        endif;

        if u_i .le stat;

            x[j]  = cand;
            cand0 = cand;

        else;

            x[j]  = cand0;
            cand0 = cand0;

        endif;

    endif;

endfor;

x;
histp(x,50);

end;