生成复制任意分布的随机数

Question

我有一组数据，其中有一个变量 z，它包含大约4000个值(范围从0.0到1.0)，这些值的直方图看起来是这样的。

现在我需要生成一个随机变量，叫做 random_z，它应该和上面的分布相似。

到目前为止，我尝试生成一个以1.0为中心的正态分布，这样我可以去掉所有大于1.0的值，从而得到一个相似的分布。我使用了 numpy.random.normal，但问题是我无法将范围限制在0.0到1.0之间，因为通常正态分布的均值是0.0，标准差是1.0。

有没有其他方法可以在Python中生成这种分布呢？

Answer 1

你可以使用拒绝采样的方法：你生成一对(z,y)，其中0<=y<=max(f(z))，一直生成直到找到一对满足y<=f(z)的值。这个生成的随机数就是z。

这个方法的好处是可以用于任何分布，但可能需要很多次尝试才能找到有效的(z,y)对。

Answer 2

如果你能大致计算出某个分布的累积分布函数（比如通过对直方图进行累加），那么从这个分布中抽样就变得非常简单了。

Sample uniformly p in interval [0.0,1.0]
Lookup the value of x at which cdf(x) == p

我想这基本上就是涉及到Pandas的那个答案在做的事情。

Answer 3

在使用 numpy.random.normal 这个函数时，你可以传入一些关键词参数来设置你得到的数组的平均值和标准差。这些关键词参数分别是 loc（表示平均值）和 scale（表示标准差）。

import numpy as np
import matplotlib.pyplot as plt

N = 4000
mean = 1.0
std = 0.5
x = []

while len(x) < N:
    y = np.random.normal(loc=mean, scale=std, size=1)[0]
    if 0.0 <= y <= 1.0:
        x.append(y)

plt.hist(x)
plt.show()

Answer 4

如果你想要进行引导采样，可以使用 random.choice() 在你观察到的序列上进行操作。

在这里，我假设你想要平滑处理，而不太关心生成新的极端值。

可以使用 pandas.Series.quantile() 和一个均匀分布的 [0,1] 随机数生成器，具体步骤如下。

训练阶段

• 把你的随机样本放入一个 pandas Series，称这个序列为 S

生产阶段

1. 生成一个随机数 u，范围在 0.0 到 1.0 之间，通常可以用 random.random() 来实现。
2. 返回 S.quantile(u)

如果你更喜欢使用 numpy 而不是 pandas，根据快速阅读的结果，第二步可以用 numpy.percentile() 来替代。

工作原理：

从样本 S 中，使用 pandas.series.quantile() 或 numpy.percentile() 来计算逆累积分布函数，这个方法叫做 逆变换采样。量化函数或百分位函数（相对于 S）将一个均匀的 [0,1] 伪随机数转换为一个具有样本 S 范围和分布的伪随机数。

简单示例代码

如果你想减少编码工作，不想写只返回单个结果的函数，那么 numpy.percentile 可能比 pandas.Series.quantile 更合适。

假设 S 是一个已有的样本。

u 将是新的均匀随机数。

newR 将是从类似 S 的分布中抽取的新随机数。

>>> import numpy as np

我需要一个随机数样本来放入 S 中。

为了创建一个示例，我将一些均匀的 [0,1] 随机数的三次方作为样本 S。通过这种方式生成示例样本，我可以提前知道——根据 (x^3)(dx) 从 0 到 1 的定积分——S 的均值应该是 1/(3+1) = 1/4 = 0.25。

在你的应用中，你可能需要做其他事情，比如读取一个文件，来创建一个包含要复制分布的数据样本的 numpy 数组 S。

>>> S = pow(np.random.random(1000),3)  # S will be 1000 samples of a power distribution

在这里，我将检查 S 的均值是否为 0.25，如上所述。

>>> S.mean()
0.25296623781420458 # OK

获取最小值和最大值，以展示 np.percentile 的工作原理。

>>> S.min()
6.1091277680105382e-10
>>> S.max()
0.99608676594692624

numpy.percentile 函数将 0-100 映射到 S 的范围。

>>> np.percentile(S,0)  # this should match the min of S
6.1091277680105382e-10 # and it does

>>> np.percentile(S,100) # this should match the max of S
0.99608676594692624 # and it does

>>> np.percentile(S,[0,100])  # this should send back an array with both min, max
[6.1091277680105382e-10, 0.99608676594692624]  # and it does

>>> np.percentile(S,np.array([0,100])) # but this doesn't.... 
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line 2803, in percentile
    if q == 0:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

如果我们生成 100 个新值，从均匀分布开始，这样做效果不太好：

>>> u = np.random.random(100)

因为这样会出错，u 的范围是 0-1，而需要的是 0-100。

这样做就可以：

>>> newR = np.percentile(S, (100*u).tolist()) 

这样可以正常工作，但如果你想要返回一个 numpy 数组，可能需要调整其类型。

>>> type(newR)
<type 'list'>

>>> newR = np.array(newR)

现在我们有了一个 numpy 数组。让我们检查一下新随机值的均值。

>>> newR.mean()
0.25549728059744525 # close enough