Welcome, guest | Sign In | My Account | Store | Cart

Generating N random numbers that probability distribution fits to any given function curve.

Python, 44 lines
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
# Generating N random numbers that probability distribution
# fits to any given function curve
# FB - 201006137

import math
import random

# define any function here!
def f(x):
    return math.sin(x)
# f(x) = 1.0 : for uniform probability distribution

# f(x) = x : for triangular probability distribution
# (math.sqrt(random.random()) would also produce triangular p.d. though.)

# f(x) = math.exp(-x*x/2.0)/math.sqrt(2.0*math.pi) : for std normal p.d.
# (taking average of (last) 2,3,... random.random() values would also
# produce normal probability distributions though.)

# define any xmin-xmax interval here! (xmin < xmax)
xmin = 0.0
xmax = math.pi

# find ymin-ymax
numSteps = 1000000 # bigger the better but slower!
ymin = f(xmin)
ymax = ymin
for i in range(numSteps):
    x = xmin + (xmax - xmin) * float(i) / numSteps
    y = f(x)
    if y < ymin: ymin = y
    if y > ymax: ymax = y

n = 10 # how many random numbers to generate
for i in range(n):
    while True:
        # generate a random number between 0 to 1
        xr = random.random()
        yr = random.random()
        x = xmin + (xmax - xmin) * xr
        y = ymin + (ymax - ymin) * yr
        if y <= f(x):
            print xr
            break

2 comments

Pietro Berkes 13 years, 10 months ago  # | flag

A better and faster way to compute random number with arbitrary distribution is to draw a number x between 0 and 1 and return cdf^{-1}(x), where cdf^{-1} is the inverse cumulative distribution function of 'f'.

FB36 (author) 13 years, 10 months ago  # | flag

The problem w/ that is how to find that function for any given arbitrary probability distribution function. I still think this way is much simpler and practical.