3

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
}
Community
  • 1
  • 1
Simon Eismann
  • 273
  • 1
  • 5
  • 17

3 Answers3

5

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
Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
  • Thank you, this works fine. Can you explain why i need to inverse the cdf? – Simon Eismann Dec 27 '16 at 10:42
  • ohhhh now i get it, x <- runif(n) ifelse(x < 0.5, 2 * x + 1, 4 * x) these are two lines (i keep forgetting i need no ; in R) You first generate a uniform [0,1] and then map it. Thank you again – Simon Eismann Dec 27 '16 at 10:50
1

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.

Community
  • 1
  • 1
Glen_b
  • 7,883
  • 2
  • 37
  • 48
  • Thank you, your approach seems very intuitive if all "bins" have the same width. I am having trouble understanding your line sample(c(1,1,2,3),n,replace=TRUE)+runif(n) – Simon Eismann Dec 27 '16 at 10:46
  • Btw you can use probabilities with the sample function like this: sample(c(1, 2, 3), size=3000000, replace=TRUE, prob=c(0.5, 0.25, 0.25)) then you dont need to use the workaround with the two 1s :-) – Simon Eismann Dec 27 '16 at 10:48
  • The (1,1, ... part is because 1 bins are twice as common; if sample is implemented well, it should be faster that way than the more general probability weighting. The ...2,3) part is splitting the 2-4 bin so they they're all the same width. Again, this is for speed. If the bin heights and widths are not all rational numbers (though nothing in the question suggests this), you will want a slightly slower but more general approach that can still be done by appropriate use of the `sample` and `runif` functions – Glen_b Dec 28 '16 at 00:17
  • Thank you for explaining, that makes sense :-) – Simon Eismann Dec 28 '16 at 14:00
0

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

Roasty247
  • 679
  • 5
  • 20