Numpy.random.normal不同分布:从分布中选择值

Question

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

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中进行采样实际上是相当困难的。关于如何有效而准确地从标准分布族中抽样，有大量的书籍。

对于你给出的例子，你可能会使用一个自定义的反转方法。

Answer 2

如果你想从任意分布中抽样，你需要累积密度函数的逆(不是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()

Answer 3

为什么不使用eval并将分布放在字符串中?

>>> cmd = "numpy.random.normal(500)"
>>> eval(cmd)

您可以根据需要操作字符串来设置分布