Generate random number from custom distribution

Question

I am trying to generate random numbers from a custom distribution, i already found this question: Simulate from an (arbitrary) continuous probability distribution but unfortunatly it does not help me since the approach suggested there requires a formula for the distribution function. My distribution is a combination of multiple uniform distributions, basically the distribution function looks like a histogram. An example would be:

f(x) = { 
    0     for  x < 1
    0.5   for  1 <= x < 2
    0.25  for  2 <= x < 4
    0     for  4 <= x
}

Answer 1

You just need inverse CDF method:

samplef <- function (n) {
  x <- runif(n)
  ifelse(x < 0.5, 2 * x + 1, 4 * x)
  }

Compute CDF yourself to verify that:

F(x) = 0                 x < 1
       0.5 * x - 0.5     1 < x < 2
       0.25 * x          2 < x < 4
       1                 x > 4

so that its inverse is:

invF(x) = 2 * x + 1      0 < x < 0.5
          4 * x          0.5 < x < 1

Answer 2

You can combine various efficient methods for sampling from discrete distributions with a continuous uniform.

That is, simulate from the integer part Y=[X] of your variable, which has a discrete distribution with probability equal to the probability of being in each interval (such as via the table method - a.k.a. the alias method), and then simply add a random uniform [0,1$, X = Y+U.

In your example, you have Y taking the values 1,2,3 with probability 0.5,0.25 and 0.25 (which is equivalent to sampling 1,1,2,3 with equal probability) and then add a random uniform.

If your "histogram" is really large this can be a very fast approach.

In R you could do a simple (if not especially efficient) version of this via

sample(c(1,1,2,3))+runif(1)

or

sample(c(1,1,2,3),n,replace=TRUE)+runif(n)

and more generally you could use the probability weights argument in sample.

If you need more speed than this gets you (and for some applications you might, especially with big histograms and really big sample sizes), you can speed up the discrete part quite a bit using approaches mentioned at the link, and programming the workhorse part of that function in a lower level language (in C, say).

That said, even just using the above code with a considerably "bigger" histogram -- dozens to hundreds of bins -- this approach seems - even on my fairly nondescript laptop - to be able to generate a million random values in well under a second, so for many applications this will be fine.

Answer 3

If the custom distribution from which you are drawing a random number is defined by empirical observations, then you can also just draw from an empirical distribution using the fishmethods::remp() package and function.

https://rdrr.io/cran/fishmethods/man/remp.html