new;

library pgraph;graphset;

format /rd 10,8;

rndseed 1166;

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

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

B = 5;

F0 = 0;

for j (1, n, 1);

    y0 = rndu(1,1)*5;

    for burn (1, 15, 1);

        F1 = 1-exp(-y0*B);
        uhat = rndu(1,1)*F1;

        x_i = -ln(1-uhat)/y0;

        F1 = 1-exp(-x_i*B);
        uhat = rndu(1,1)*F1;

        y_1 = -ln(1-uhat)/x_i;

        y0 = y_1;

    endfor;

    x[j] = x_i;

endfor;

x;
histp(x,50);

end;