CODEWITHSUNDEEP
pythonrandomscipydistribution
rfeynman
rfeynman

Reputation: 121

Generate Random variates with a certain Probability density function by python

I want to generate a random variate with maxwell energy(e.g. f(ene)) or arbitrary distribution. The probablily density function is like :

def f(ene):  
     le=3  
     return 2(ene/pi/le**3)*np.exp(-ene/le)

I would like to generate 10000 samples like

f.rvs(scale,size)

which generate a list includes 10000 elements which density function is match on f(ene)

How could I do it?

Example: http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.maxwell.html scipy provide a maxwell distribution: with method: rvs(loc=0, scale=1, size=1) Random variates.

print maxwell.rvs(1,10000)

will generate 10000 sample with maxwell distribution. But scipy only provide some functions. Now I have another function not in scipy's . How could I do that?

Upvotes: 0

Views: 2287

Answers (3)

Christoph Baumgarten
Christoph Baumgarten

Reputation: 189

In the past few years, nice new tools have been added to SciPy to address this kind of problem in Python. You can easily generate samples from custom continuous or discrete univariate distributions by just providing some information about the distribution, such as the density / pdf.

Here is an example:

import numpy as np
from scipy.stats.sampling import NumericalInversePolynomial
from matplotlib import pyplot as plt
from scipy.integrate import quad


class MyDist:
    def __init__(self, a):
        self.a = a

    def support(self):
        # distribution restricted to 0, 5, can be changed as needed
        return (0, 5)

    def pdf(self, x):
        # this is not a proper pdf, the normalizing
        # constant is missing (does not integrate to one)
        return x * (x + np.sin(5*x) + 2) * np.exp(-x**self.a)


dist = MyDist(0.5)
gen = NumericalInversePolynomial(dist)

# compute the missing normalizing constant to plot the pdf
const_pdf = quad(dist.pdf, *dist.support())[0]

r = gen.rvs(size=50000)
x = np.linspace(r.min(), r.max(), 500)

# show histogram together with the pdf
plt.plot(x, dist.pdf(x) / const_pdf)
plt.hist(r, density=True, bins=100)
plt.show()

Samples and density

Overview of the different methods: https://docs.scipy.org/doc/scipy/reference/stats.sampling.html

Tutorial: https://docs.scipy.org/doc/scipy/tutorial/stats/sampling.html

Upvotes: 1

Severin Pappadeux
Severin Pappadeux

Reputation: 20080

Ok, lets start with distribution.

Probability density function for speed is

PDF(v) = sqrt(2/PI)*v^2*exp(-v^2/2a^2)/a^3

Lets convert it from speed to energy, considering E=v^2/2 (and assume mass is equal to 1)

PDF(E) = sqrt(2/PI)*2*E*exp(-E/a^2)/a^3

Looks ok? No, it is wrong. Why? Because PDF times the interval gives you probability, and we didn't convert interval. Probability for speed was PDF(v) * dv, and we have to insert Jacobian

PDF(E)*dE = sqrt(2/PI)*2*E*exp(-E/a^2)/a^3 * |dv/dE| * dE

|dv/dE| = 1/|dE/dv| = 1/v = 1 / sqrt(2*E)

All together

PDF(E) = 2*sqrt(1/PI)*sqrt(E)*exp(-E/a^2)/a^3

where a=sqrt(k*T)

Your distribution is missing sqrt(E)

How I would sample it. I would use the fact that it is actually product of three independent 1D distribution

PDF(vx) = 1/sqrt(2*PI)*exp(-vx^2/2a^2)
PDF(vy) = 1/sqrt(2*PI)*exp(-vy^2/2a^2)
PDF(vz) = 1/sqrt(2*PI)*exp(-vz^2/2a^2)

which means they are all Gaussian. So simple algorithm: sample vx from gaussian with mean equal to 0 and sigma equal to a, sample vy independently from the same distribution, sample vz and combine them all together

v = sqrt(vx^2 + vy^2 + vz^2)

or if you need sampled energy

E = ( vx^2 + vy^2 + vz^2 )/ 2

Upvotes: 1

James Mertz
James Mertz

Reputation: 8759

This is the best that I can come up with, but I'm not confident on what you're trying to do here:

import numpy as np

def f(ene):
    le = 3
    return 2 * (ene / np.pi / le ** 3) * np.exp(-ene / le)

def rvs(scale, size):
    samples = np.random.sample((size,))
    return scale * f(samples)

I got this as a result. Obviously this will be different each time as this is randomized. 

In:  print(rvs(100, 10))
Out: [ 1.56035154  0.85509302  0.76543496  1.36966862  0.20596924  0.0395071
  1.35029318  0.69437599  0.59772671  0.56721737]

Upvotes: 0

Related Questions