How to distribute values to have a specified probability density function

Question

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?

I typed "random sampling with distribution for fixed sum" in google, and got a ton of results, especially questions asked at https://stats.stackexchange.com ; for instance, https://stats.stackexchange.com/questions/315109/draw-from-random-distributions-so-that-sum-of-values-is-fixed https://stats.stackexchange.com/questions/220543/generate-a-random-set-of-numbers-with-fixed-sum-and-desired-means-and-variances https://stats.stackexchange.com/questions/243260/sample-random-variables-conditional-on-their-sum — Stef, Sep 20 '21 at 13:33
https://stats.stackexchange.com/questions/23955/how-to-sample-natural-numbers-such-that-the-sum-is-equal-to-a-constant https://stackoverflow.com/questions/5622608/choosing-n-numbers-with-fixed-sum https://stackoverflow.com/questions/29187044/generate-n-random-numbers-within-a-range-with-a-constant-sum https://stackoverflow.com/questions/53279807/how-to-get-random-number-list-with-fixed-sum-and-size https://stackoverflow.com/questions/49016047/random-sampling-to-give-an-exact-sum — Stef, Sep 20 '21 at 13:35
This question in particular appears to be exactly what you are asking: [Sample random variables conditional on their sum](https://stats.stackexchange.com/questions/243260/sample-random-variables-conditional-on-their-sum) — Stef, Sep 20 '21 at 13:47
@JoBi Please note that I found an error and fixed it before you use my code. — Lukas S, Sep 21 '21 at 10:59

score 3 · Answer 1 · edited Sep 20 '21 at 16:49

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)

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:

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

Lukas S · Accepted Answer · 2021-09-21T10:57:00.217

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

How to distribute values to have a specified probability density function

2 Answers2