How to sample from a given probability distribution?

Question

I am plotting a histogram with the same probability distribution of particle in a ground state one dimensional box, using random numbers. But when compared with the original distribution, the values are getting cut at the top.

Here is the code

from numpy import*
from matplotlib.pyplot import*

lit=[]
def f(x):
    return 2*(sin(pi*x)**2)
for i in range(0,10000):
    x=random.uniform(0,1)
    y=random.uniform(0,1)
    if y<f(x):
        lit.append(x)

l=linspace(0,1,10000)
hist(lit,bins=100,density=True)
plot(l,f(l))
show()

The graph produced is:

How to improve this code to match the original?

Does this answer your question? [How to generate random numbers with predefined probability distribution?](https://stackoverflow.com/questions/51050658/how-to-generate-random-numbers-with-predefined-probability-distribution) — Peter O., Jan 18 '21 at 10:44
See also: https://stackoverflow.com/questions/55724501/generate-a-random-number-according-to-a-specific-distribution-with-a-function/64961053#64961053 — Peter O., Jan 18 '21 at 10:45

Julien · Answer 1 · 2021-01-18T07:57:27.030

Your issue is that your density function f ranges in [0,2] but you draw y from [0,1].

Using y=np.random.uniform(0,2) will fix it.

But a better approach is to remap the uniform density to your desired function directly:

from scipy.interpolate import interp1d

l = np.linspace(0,1,100)
g = interp1d(l-np.sin(2*np.pi*l)/2/np.pi, l) # inverse for integral of f
# (not sure if analytical solution to above exists... so using interpolation instead)
lit = g(np.random.rand(10000))

In the general case where f may not have an analytical integral, you can use numerical integral too:

dl = l[1]
g = interp1d(np.cumsum(f(l))*dl, l)

How to sample from a given probability distribution?

1 Answers1