Algorithm:
1. Set m=0
2. Determine Py_max and the range (ymin and ymax) of the target distribution, Py(y).
3. Generate N random numbers, x(N), y(N), from a uniform distribution.
4. For each value of y: Set y = ymin + (ymax-ymin)*y
5. LOOP from i=1 to N
Generate a new random number, ptest distributed
uniformly in the range 0 to Py_max:
Set ptest = Py_max * y(i)
Test: if(Py( y(i) ) > ptest) then
Set m=m+1 # Entry y(i) passed the test
Set z(m)=y(i) # so we record it in z(m)
otherwise
remove(reject) y(i) from the squence
end if Test
END LOOP i
5. The numbers remaining in the sequence are therefore distributed according
to the function Py(y).
Plot z
|
