CODEWITHSUNDEEP
rrandomstatisticsprobability
Chris
Chris

Reputation: 41

Generating a random sample for a given probability distribution

I have a random variable X with pdf f(x)=4xe^-x, where x>0.

How do i draw a random sample of size, let's say, 1000 from this distribution?

Upvotes: 0

Views: 1598

Answers (2)

duckmayr
duckmayr

Reputation: 16910

One easy solution is rejection sampling (though see my comments on Severin Pappadeux's answer below). A general rejection sampling algorithm is pretty easy to implement in R (and I would be surprised if it were not implemented in several packages, I just didn't look up which packages would have such functionality):

#' Rejection sampling
#' 
#' @param f A function calculating the density of interest
#' @param g A function giving the proposal density
#' @param rg A function providing random samples from g
#' @param M A numeric vector of length one giving a bound on the the ratio
#'   f(x) / g(x). M must be > 1 and greater than or equal to f(x) / g(x) over
#'   the whole support of X.
#' @param n An integer vector of length one giving the number of samples to
#'   draw; the default is ten thousand.
#' @param ... Further arguments to be passed to g and rg
rejection_sampling <- function(f, g, rg, M, n = 10000, ...) {
    result <- numeric(n)
    for ( i in 1:n ) {
        reject <- TRUE
        while ( reject ) {
            y <- rg(n = 1, ...)
            u <- runif(1)
            if ( u < ( f(y) / (M * g(y, ...)) ) ) {
                result[i] <- y
                reject <- FALSE
            }
        }
    }
    return(result)
}

Then we can use that to get a sample from your distribution, and plot the sample's density along with the true probability density to see how well it worked:

x <- seq(0.01, 15, 0.01)
f <- function(x) 4 * x * exp(-x)
y <- f(x)
set.seed(123)
z <- rejection_sampling(f = f, g = dexp, rg = rexp, M = 10, n = 1e3, rate = 1/4)
dens <- density(z, from = 0.01, to = 15)
scaling_constant <- max(y) / max(dens$y)
plot(x, y, type = "l", xlab = "x", ylab = "f(x)", lty = 2, col = "blue")
lines(dens$x, dens$y * scaling_constant, col = "red", lty = 3)
legend("topright", bty = "n", lty = 2:3, col = c("blue", "red"),
       legend = c("True f(x)", "(Re-scaled) density of sample"))

enter image description here

Rejection sampling works by taking samples from a proposal distribution and rejecting them if a random uniform deviate is greater than the ratio f(x) / M g(x), where g(x) is your proposal density, and M is a bound on f(x) / g(x) as briefly described in the Roxygen documentation above.

I used above an exponential proposal distribution with a rate parameter of 1/4. You could use others.

This p.d.f. is proportional to a Gamma distribution with a shape of 2 and scale of 1, as mentioned by Severin Pappadeux's answer. (That is, if you plug those parameters into the Gamma p.d.f., you'll see it differs from yours only by a scaling constant). Depending on what you're doing this for, that could be a better way to go, or this way. I'm not sure if your objective is to generate samples from arbitrary distributions such as this one as an example, or you need samples from this distribution itself, etc.... Usually it's best to recognize if your arbitrary p.d.f. is actually an example of an implemented distribution, but if you're not in that situation, you need something like rejection sampling.

Upvotes: 3

Severin Pappadeux
Severin Pappadeux

Reputation: 20080

This is Gamma distribution with shape=a=2 and scale=1.

q <- rgamma(1000, shape = 2, scale = 1)

Upvotes: 4

Related Questions