There is theory for rejection sampling that tells you how long you will
need to wait: see any good book on *simulation*.
However,
1) Your example runs in a moderate number of simulations for me, but for
some argument values it can be arbitrarily slow.
2) There are much better envelopes, for example a normal with a shifted
mean, or a t.
On Sun, 25 May 2003, Stefano Sofia wrote:
> Dear Splus user, let me share with you a sampling problem.
>
> In order to sampling from an unknown distribution g,
That's impossible! You have to know the pdf up to a constant to use
rejection sampling.
> I'm trying to implement
> the Rejection Sampling procedure explained by Gilks, Richardson and
> Spiegelhalter in the book
> "Markov Chain Monte Carlo in Practice".
> The key points of this procedure are two:
> 1. an envelope function G is required, where for every value Y1 sampled from
> G, G(Y1) > g(Y1)
> 2. sampled a value Y1 from G, Y1 can be accepted only if, sampling a value U
> from a Uniform(0,1), U <= g(Y1)/G(Y1)
You can't do that. If g and G are densities (and they should be), it is
impossible to have two pdfs with g < G everywhere. And if you did not
mean densities, what did you mean?
> In my case,
> g = f(x)P(x)
That is not a density!
> where
> 1. f(x) is a Normal distribution with a certain mean and standard deviation
> 2. P(x) = exp(Y1)/(1 + exp(Y1))
> Because of P(x) is always in (0,1) as envelope function G I can use f(x);
> nothing easier, being f(x) a Normal distribution.
To sample from pdf g, what you need is an envelope pdf G such that g <=
MG: you then accept with probability (g/MG)(Y1).
You have g' = fP \propto g, say cg' = g. Now g'/G < 1 and so
g/G < c gives a rejection method with c=M and hence accepting if
g/G > cU or g'/G = P > U. Such a method requires on average 1/M steps.
[...]
--
Brian D. Ripley, ripley@stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
|