Generating a random number based on a distribution function

Question

We can generate random numbers in an interval [a,b] easily if we want to make it uniformely:

A=rand()*(b-a) + a

where rand() is a function which can generate a uniform random number between 0 and 1. so A is a random number in [a,b].

For generating a random number based on a distribution function like y=x-x^2, i faced a problem.

I would like to use a method mentioned here. But i am not interested to use the python function inverse_cdf(np.random.uniform()).

I can compute the CDF of function "y" by an integration over 0 and X and i call it "f". But when i put the rand() function(a number between 0 and 1) into the inverse function of f, i get a complex number! It means: A=f^(-1) (rand()) returns a complex number.

Is it a correct way for generating a random number based on a distribution function?

I used this website to compute the inverse of f=x^2/2 - x^3/3 and the code below is a part of the calculation that shows that the tmp1 is always negative

for i=1:10
rnd1=rand;
tmp1  = 2*sqrt(6)*sqrt(6*rnd1^2-rnd1)-12*rnd1+1
cTmp1 = tmp1^(1/3)
end

I'm struggling to understand the question, could you add some example code and tag with the language you're using? — Sam Mason, Jan 22 '21 at 13:31
I edited the topic a bit. Actually i am going to implement it in matlab and i don't want to use any existing function. I think that my method has some problem. Maybe i misunderstood the method for generating the random number. — Denis, Jan 22 '21 at 13:44
Dunno, I get a much larger equation: https://www4c.wolframalpha.com/Calculate/MSP/MSP309811adbcfd3e7473h400006acbbaia94eidce9?MSPStoreType=image/gif&s=48&w=472.&h= — Ander Biguri, Jan 22 '21 at 15:00
Your fundamental problem is that `x - x^2` is not a valid density function. To be kosher, a density `f(x)` must satisfy two properties: 1) `f(x) ≥ 0` for all x, and 2) the integral from -infinity to infinity must be 1. You need to specify a range, scaling, and possibly an additive constant to turn your equation into a legitimate density. That's an essential starting point to deriving an inverse transform generating algorithm. — pjs, Jan 22 '21 at 15:31
To get concrete, `x - x^2` is only positive for `0 ≤ x ≤ 1`. Integrating over that range yields an area of 1/6. So you could make it a valid density by specifying `f(x) = 6(x - x^2) for 0 ≤ x ≤ 1; 0 otherwise`. But you could also add a constant to shift the equation up or down, which would affect both the range and the resulting area. So what is the actual density you want? — pjs, Jan 22 '21 at 15:36
By the way, I'd recommend that you stick with standard notation. Density functions are denoted `f`, CDFs are denoted `F`. They are two related but different functions, and using standard notation helps you communicate and keep your head straight regarding the steps involved in solving the problem. — pjs, Jan 22 '21 at 15:46
@pjs thank you for the hints. I completely forgot about the properties... — Denis, Jan 23 '21 at 12:44
@pjs I calculated the integral of the valid `f(x)` from 0 to z, and i got the CDF as `F(z)=3z^2-2z^3`. So the problem is that inverting this function, again gives me `sqrt(X^2-X)`. So if i put a random variable between 0 and 1, i get a complex number. — Denis, Jan 23 '21 at 13:41
That's the same F I came up with under the assumption that there was no additive constant and the range was [0,1]. I think you may have done your algebra incorrectly, because if you plot the CDF there is clearly a 1-1 mapping between F and x, so real-valued solutions must exist. — pjs, Jan 23 '21 at 16:04
I inverted `F(z)` with this [website](https://www.wolframalpha.com/widgets/view.jsp?id=d08726019e4a2a15cb1d49092e4d0522). There is `sqrt(Z^2-Z)` in the solution. — Denis, Jan 23 '21 at 16:41
[Wolfram doesn't give a closed form solution when I add the constraints](https://www.wolframalpha.com/input/?i=solve+3x%5E2+-+2x%5E3+%3D+U+for+0+%3C+x+%3C+1+where+0+%3C+U+%3C+1), so it resorts to plotting. You should be able to solve it numerically using Newton's method given that even Wolfram doesn't yield a clean algebraic solution. — pjs, Jan 23 '21 at 16:53

score 3 · Accepted Answer · answered Jan 22 '21 at 14:35

The question is: "Is it a correct way for generating a random number based on a distribution function?" Thus, basically you have either continuous Probability Density Function (PDF) or discrete Probability Mass Function (PMF) with notation f(x) and you are are looking for a way to find random variate x.

There are at least two ways to do this.

To use Inverse Transform Distribution.
To use Rejection method

Using Inverse transform: If we know the function of probability distribution, then for some Cumulative Distribution Function (CDF) we can find the closed from of random variate. Suppose your probability function is f(x) and the CDF is F(x) then assuming you can get the inverse function, you can get random variate

x=inverse F(U)

where U is random uniform distribution

Using Rejection Method: If the CDF does not have a closed form inverse, then you can always use rejection method. The idea of rejection method is to generate a random 2D point (a pair of random number): (r1, r2) and then the point is either under the curve or over the curve of the PDF. If the point is under the curve, then we take this value as our random number, otherwise we sample another point. Suppose PDF f(x) is bounded by a maximum value M

r1 is generated within interval [a, b] of horizontal axis
r2 is generated within interval [0, M] of vertical axis
If r2 <f (r1) then accept r1 and output r1
   Else reject r1 and generate new random point

Inverse Transform method is superior to the Rejection method if you can find the inverse of the CDF because you can get the closed form. For instance, CDF of Exponential distribution with rate 1 is F(x) = 1-exp(-x). The inverse transform would be

x= inverse F(U) = -ln(1-U) = -ln(U)

because 1-U2=U2

Generating a random number based on a distribution function

1 Answers1