numpy.random.正常不同分布：从分布中选择值

Question

我有一个能量的幂律分布，我想根据这个分布选取n个随机的能量。我尝试过用随机数手动执行此操作，但这对于我想做的事情来说效率太低了。我想知道在numpy（或其他）中是否有一个方法可以像numpy.random.normal一样工作，除了可以指定分布而不是使用正态分布。所以在我看来，一个例子可能看起来像（类似于numpy.random.正常)公司名称：

import numpy as np

# Energies from within which I want values drawn
eMin = 50.
eMax = 2500.

# Amount of energies to be drawn
n = 10000

photons = []

for i in range(n):

    # Method that I just made up which would work like random.normal,
    # i.e. return an energy on the distribution based on its probability,
    # but take a distribution other than a normal distribution
    photons.append(np.random.distro(eMin, eMax, lambda e: e**(-1.)))

print(photons)

打印photons应该会给我一个长度为10000的列表，由这个分布中的能量填充。如果我把这个柱状图，它在低能量下会有更大的bin值。在

我不确定这种方法是否存在，但它似乎应该存在。我希望我想做的事情很清楚。在

编辑：

我见过numpy.random.power，但我的指数是-1，所以我认为这行不通。在

Answer 1

如果你想从任意分布中取样，你需要累积密度函数的倒数（而不是pdf）。在

然后从[0,1]范围内均匀地抽取一个概率，并将其输入cdf的逆函数，以得到相应的值。在

通常不可能从pdf分析中获得cdf。 但是，如果你愿意近似分布，你可以通过在它的域上以固定的间隔计算f（x），然后在这个向量上做一个求和得到cdf的近似值，然后从这个近似值得到逆。在

粗略代码片段：

import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate

def f(x):
   """
   substitute this function with your arbitrary distribution
   must be positive over domain
   """
   return 1/float(x)


#you should vary inputVals to cover the domain of f (for better accurracy you can
#be clever about spacing of values as well). Here i space them logarithmically
#up to 1 then at regular intervals but you could definitely do better
inputVals = np.hstack([1.**np.arange(-1000000,0,100),range(1,10000)])

#everything else should just work
funcVals = np.array([f(x) for x in inputVals])
cdf = np.zeros(len(funcVals))
diff = np.diff(funcVals)
for i in xrange(1,len(funcVals)):
   cdf[i] = cdf[i-1]+funcVals[i-1]*diff[i-1]
cdf /= cdf[-1]

#you could also improve the approximation by choosing appropriate interpolator
inverseCdf = scipy.interpolate.interp1d(cdf,inputVals)

#grab 10k samples from distribution
samples = [inverseCdf(x) for x in np.random.uniform(0,1,size = 100000)]

plt.hist(samples,bins=500)
plt.show()