Dear Splus user, let me share with you a sampling problem.
In order to sampling from an unknown distribution g, 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)
In my case,
g = f(x)P(x)
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.
When
1. the mean of f(x) is negative
2. the standard deviation very small
then
1. Y1 is a value generally close to the mean
2. G(Y1) is a quite high value
3. g(Y1) is a quite small value
4. g(Y1)/G(Y1) = f(Y1)P(Y1)/f(Y1) = P(Y1) is a small value, close to 0
and it becomes difficult to sample a value U from a Uniform(0,1) smaller that
P(Y1).
In my calculations, values approximatively are
1. mean for the normal distribution = -5.1
2. variance for the normal distribution = 0.00002
3. g(Y1) / G(Y1) = P(Y1) = 0.005
without finding a value U that validates the sampled value Y1, not even after
8 hours of calculations.
But if I sample for example 1000 random values from a Uniform(0,1) using
runif(1000, min = 0, max = 1)
I can easily get some values smaller than 0.005.
Then, why my loop is "endless"? My pc is a P IV.
Is it due to syntax problems? (my procedure won't be te best one and not the
most optimised, but I believe it is correct)
Is it due to probability matters and therefore nothing canbe done?
Is there something I am missing or I can do?
thank you for your usual attention
hints and suggestions will be appreciated
Stefano
The original algorithm is:
REPEAT
SAMPLE A POINT Y1 FROM G
SAMPLE A UNIFORM(0,1) RANDOM VARIABLE U
IF (U <= g(Y1)/G(Y1)) {ACCEPT Y1}
UNTIL ONE Y1 IS ACCEPTED
My Splus procedure is
rejectionsampling <- function(normalmean, normalvariance)
{
repeat {
# REPEAT THIS LOOP UNTIL ONE Y IS ACCEPTED
# STEP 1: SAMPLE A POINT Y1 FROM G
Y1 <- rnorm(1, mean = normalmean, sd = sqrt(normalvariance))
# STEP 2: SAMPLE A UNIFORM(0,1) RANDOM VARIABLE U
U <- runif(1, min = 0, max = 1)
# STEP 3: IF U <= g(Y1)/G(Y1) ACCEPT Y1;
P <- exp(Y1)/(1 + exp(Y1))
if(U <= P) {
sampledvalue <- Y1
break
}
}
result <- sampledvalue
result
}
|