CODEWITHSUNDEEP
pythonrandomsampling
Ihsan Ahmed K
Ihsan Ahmed K

Reputation: 159

How to sample from a given probability distribution?

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: enter image description here

How to improve this code to match the original?

Upvotes: 1

Views: 2133

Answers (1)

Julien
Julien

Reputation: 15206

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)

Upvotes: 1

Related Questions