如何分配值以具有指定的概率密度函数
How to distribute values to have a specified probability density function
我有一个概率密度函数(pdf)如下
: f(m) = 4 m exp(-2m)
公式的图表附在下面
如何取N个粒子并在这些粒子之间分配能量N(平均能量为1),使得分布的pdf遵循上面的图形和等式?
我找到了解决这个问题的简单方法,它对我的研究非常有效。我使用 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
。
编辑:在我的计算中发现了一个讨厌的错误但已修复。
我将使用 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()
我有一个概率密度函数(pdf)如下
公式的图表附在下面
如何取N个粒子并在这些粒子之间分配能量N(平均能量为1),使得分布的pdf遵循上面的图形和等式?
我找到了解决这个问题的简单方法,它对我的研究非常有效。我使用 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
。
编辑:在我的计算中发现了一个讨厌的错误但已修复。
我将使用 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()