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的列表，里面填满了这个分布下的能量值。如果我把这些值做成直方图，低能量的区间会有更多的值。

我不太确定是否有这样的函数，但感觉应该有。我希望我说的能让人明白我想做什么。

补充：

我看到过numpy.random.power，但我的指数是-1，所以我觉得这个方法不适用。

Answer 1

你为什么不使用 eval 呢，把分布放在一个字符串里？

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

你可以随意处理这个字符串来设置分布。

Answer 2

如果你想从一个任意的分布中抽样，你需要用到累积分布函数的反函数（而不是概率密度函数）。

接着，你可以在[0,1]这个范围内均匀地抽取一个概率值，然后把这个值放入累积分布函数的反函数中，就能得到对应的结果。

通常情况下，从概率密度函数得到累积分布函数并不容易。不过，如果你愿意对分布进行近似，你可以在分布的范围内定期计算f(x)的值，然后对这些值进行累加，得到一个累积分布函数的近似值，接着再从这个近似值中求出反函数。

下面是一个简单的代码示例：

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

从任意的概率分布中进行抽样其实是挺难的。关于如何高效且准确地从标准分布中抽样，有很多厚厚的书籍专门讨论这个问题。

看起来，你可以尝试用一种自定义的反转方法来处理你给出的例子。