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 中很好地采样实际上非常困难。 large and dense books 只是关于如何从标准分布族中高效准确地采样。

对于您提供的示例，您可能可以使用自定义反转方法。

Answer 2

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

然后，您从 [0,1] 范围内均匀地采样一个概率，并将其输入 cdf 的倒数，以获得相应的值。

通常无法通过解析方式从 pdf 中获取 cdf。 但是，如果您愿意近似分布，您可以通过在其域上定期计算 f(x) 来实现，然后对该向量进行 cumsum 以获得 cdf 的近似值，并由此近似得到逆。

粗略代码sn-p：

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)

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