CODEWITHSUNDEEP
numpymatplotlibmathprobability-density
Jo Bi
Jo Bi

Reputation: 97

How to distribute values to have a specified probability density function

I have a probability density function (pdf) as follows

enter image description here: f(m) = 4 m exp(-2m)

The graph for the equation is attached below

enter image description here

How can I take N particles and distribute energy N (the average energy is 1) among these particles such that the pdf of the distribution follows the above graph and equation?

Upvotes: 1

Views: 286

Answers (2)

Lukas S
Lukas S

Reputation: 3583

Edit: found a pesky bug in my calculations but fixed it.

I am going to do this using scipy.stats.rvs_ratio_uniforms. It requires me to calculate some values first. See the documentation I have linked above. We need:

umax = sup sqrt(pdf(x))
vmin = inf (x - c) sqrt(pdf(x))
vmax = sup (x - c) sqrt(pdf(x))

I will calculate them using sympy.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import rvs_ratio_uniforms
from sympy import *
import sympy 

x = symbols('x')
c = 0 # included since it could also be chosen defferently

pdf_exp = 4* x* sympy.exp(-2*x)
pdf = lambdify(x,pdf_exp,'numpy')

umax = float(maximum(sqrt(pdf_exp), x))
vmin = float(minimum((x - c) * sqrt(pdf_exp), x))
vmax = float(maximum((x - c) * sqrt(pdf_exp), x))

To make sure everything went well I generate 10**6 random samples and plot a histogram against the pdf

data = rvs_ratio_uniforms(pdf, umax, vmin, vmax, size=10**6, c=c)

t = np.linspace(0,10,10**5)
_ = plt.hist(data, bins='auto', density=True)
plt.plot(t, pdf(t))
plt.show()

histogram to check the result

Upvotes: 3

Jo Bi
Jo Bi

Reputation: 97

I found a simple solution to this question and it worked really well for my research. I made use of the numpy.random.choice function and I inserted the pdf into the probability part of the code. However, I had to change the pdf such that the sum = 1.

def initial(num, pr, xl):
    x=[]
    
    for i in range(num):
        
        x.append(np.random.choice(xl, p = pr ))
        
    return x


n = 2
a = 4
f = []
xl = np.linspace(0,5,1000)
for i in xl:
        fun = a*(i**(n-1))*np.exp(-n*i) #calculating the pdf
        f.append(fun)
        
f = [i/sum(f) for i in f] #to make sure that the sum of the probabilities is 1. 

xl = initial(10**6, np.array(f), xl) #xl is the required distribution

x1,y1 = np.histogram(xl, 50, density = True)
y1 = y1[:-1]
#to check the result
plt.plot(y1,x1)

enter image description here

One thing to note here is that I have not done anything in the code to make sure the total energy is constant. The reason for this is, in the equation

Equation:

enter image description here

Since the average energy in this case is 1, the total energy will automatically be N.

Upvotes: 3

Related Questions