如何分配值以具有指定的概率密度函数答案

【问题标题】：How to distribute values to have a specified probability density function如何分配值以具有指定的概率密度函数
【发布时间】：2021-11-14 04:14:36
【问题描述】：

我有一个概率密度函数（pdf）如下

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

方程的图表附在下面

如何获取 N 个粒子并在这些粒子之间分配能量 N（平均能量为 1），使得分布的 pdf 遵循上图和方程？

【问题讨论】：

我在 google 中输入了“随机抽样与固定总和分布”，得到了大量结果，尤其是在 stats.stackexchange.com 提出的问题；例如，stats.stackexchange.com/questions/315109/…stats.stackexchange.com/questions/220543/…stats.stackexchange.com/questions/243260/…
stats.stackexchange.com/questions/23955/… stackoverflow.com/questions/5622608/… stackoverflow.com/questions/29187044/… stackoverflow.com/questions/53279807/… stackoverflow.com/questions/49016047/…
这个问题似乎正是您要问的：Sample random variables conditional on their sum
@JoBi 请注意，在您使用我的代码之前，我发现了一个错误并修复了它。

标签： numpy matplotlib math probability-density

【解决方案1】：

我为这个问题找到了一个简单的解决方案，它对我的研究非常有效。我使用了numpy.random.choice 函数，并将pdf 插入代码的概率部分。但是，我不得不更改 pdf 以使 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)

这里要注意的一点是，我没有在代码中做任何事情来确保总能量是恒定的。原因是，在等式中

方程：

由于本例中的平均能量为1，因此总能量将自动为N。

【讨论】：

【解决方案2】：

编辑：在我的计算中发现了一个讨厌的错误，但已修复。

我将使用scipy.stats.rvs_ratio_uniforms 来执行此操作。它需要我先计算一些值。请参阅我上面链接的文档。我们需要：

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

我将使用 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))

为了确保一切顺利，我生成了10**6 随机样本并针对 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()

【讨论】：