[R] Question about acception rejection sampling - NOT R question

[Charles C. Berry]

On Fri, 13 Jul 2007, Leeds, Mark (IED) wrote:

> This is not related to R but I was hoping that someone could help me. I
> am reading the "Understanding the Metropolis Hastings Algorithm"
> paper from the American Statistician by Chip and Greenberg, 1995, Vol
> 49, No 4. Right at the beginning they explain the algorithm for basic
> acceptance  rejection sampling in which you want to simulate a density
> from f(x) but it's not easy and you are able
> to generate from another density called h(x). The assumption is that
> there exists some c such that f(x) <= c(h(x) for all x
>
> They clearly explain the steps as follows :
>
> 1) generate Z from h(x).
>
> 2) generate u from a Uniform(0,1)
>
> 3) if u is less than or equal to f(Z)/c(h(Z) then return Z as the
> generated RV; otherwise reject it and try again.
>
> I think that, since f(Z)/c(h(z) is U(0,1), then u has the distrbution as
> f(Z)/c(h(Z).
>
> But, I don't understand why the generated and accepted Z's have the same
> density  as f(x) ?
>
> Does someone know where there is a proof of this or if it's reasonably
> to explain, please feel free to explain it.

The original reference to J. von Neumann's work, which no doubt has a 
proof, is in

 	http://en.wikipedia.org/wiki/Rejection_sampling

along with a suggestive graph.

Here is another graph that may give your intuition a boost:

f <- function(x) dnorm(x)/2+dnorm(x,sd=0.1)/2
h <- dnorm
Z <- rnorm(1000000)
u <- runif(1000000)
cee <- f(0)/h(0) # as is obvious
u.lt.f.over.ch <- u < f(Z)/cee/h(Z)
cutpts <- seq(-5,5,by=0.1)
midpts <- head(cutpts,n=-1)/2 + tail(cutpts,n=-1)/2
tab <- table( !u.lt.f.over.ch, cut( Z, cutpts ) )
bp <- barplot( tab )
lines( bp, f(midpts)*sum(u.lt.f.over.ch)*0.1,col='red',lwd=2)

Of course, there is some discreteness in this plot.
If you have a 1280x1024 or wider screen or want to zoom in on a pdf(), 
you might try 0.025 or even 0.0125 in place of 0.1.

Turn this graph on its side and think about what u.lt.f.over.ch is doing 
conditionally on Z, i.e. look at one bar and think about why it is split 
where it is.


HTH,

Chuck

> They authors definitely believe it's too trivial because they don't. The
> reason I ask is because, if I don't understand this then
> I definitely  won't understand the rest of the paper because it gets
> much more complicated.  I willing to track down the proof but I don't
> know where to look. Thanks.
> --------------------------------------------------------
>
> This is not an offer (or solicitation of an offer) to buy/se...{{dropped}}
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901