[FIXED] How to distribute values to have a specified probability density function

Issue

I have a probability density function (pdf) as follows

: f(m) = 4 m exp(-2m)

The graph for the equation is attached below

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?

Solution

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()

Answered By - user2640045