从概率分布生成随机变量

Question

我已经从我的 python 数据集中提取了一些变量，我想从我的分布中生成一个更大的数据集。问题是我试图在保持类似行为的同时向新数据集引入一些可变性。这是我提取的数据的一个示例，其中包含 400 个观察值：

Value    Observation Count     Ratio of Entries
1        352                    0.88
2        28                     0.07
3        8                      0.02
4        4                      0.01
7        4                      0.01
13       4                      0.01

现在，我正尝试使用此信息生成具有 2,000 个观测值的类似数据集。我知道 numpy.random.choice 和 random.choice 函数，但我不想使用完全相同的分布。相反，我想根据分布生成随机变量（值列），但具有更多可变性。我希望更大的数据集看起来像这样的示例：

Value         Observation Count        Ratio of Entries
1             1763                     0.8815
2             151                      0.0755
3             32                       0.0160
4             19                       0.0095
5             10                       0.0050
6             8                        0.0040
7             2                        0.0010
8             4                        0.0020
9             2                        0.0010
10            3                        0.0015
11            1                        0.0005
12            1                        0.0005
13            1                        0.0005
14            2                        0.0010
15            1                        0.0005

因此，如果我用指数衰减函数拟合原始数据，就可以估计新分布，但是，我对连续变量不感兴趣。我该如何解决这个问题，是否有与我正在尝试做的事情相关的特定或数学方法？

Answer 1

听起来你想根据第二个table中描述的PDF生成数据。 PDF 类似于

0 for x <= B
A*exp(-A*(x-B)) for x > B

A 定义分布的宽度，它始终被归一化为面积 1。B 是水平偏移，在您的情况下为零。您可以通过 ceil.

装箱使其成为整数分布

归一化衰减指数的 CDF 是 1 - exp(-A*(x-B))。通常，自定义分布的一种简单方法是生成统一数字并通过 CDF 对其进行映射。

幸运的是，您不必这样做，因为 scipy.stats.expon already provides the implementation you are looking for. All you have to do is fit to the data in your last column to get A (B is clearly zero). You can easily do this with curve_fit。请记住 A 在 scipy PDF 语言中映射到 1.0/scale。

这是一些示例代码。我通过计算整数输入的 objective 函数从 n-1 到 n 的积分，在此处增加了额外的复杂性，并在进行拟合时考虑了装箱。

import numpy as np
from scipy.optimize import curve_fit
from scipy.stats import expon

def model(x, a):
    return np.exp(-a * (x - 1)) - exp(-a * x)
    #Alternnative:
    # return -np.diff(np.exp(-a * np.concatenate(([x[0] - 1], x))))

x = np.arange(1, 16)
p = np.array([0.8815, 0.0755, ..., 0.0010, 0.0005])
a = curve_fit(model, x, p, 0.01)
samples = np.ceil(expon.rvs(scale=1/a, size=2000)).astype(int)
samples[samples == 0] = 1
data = np.bincount(samples)[1:]

Answer 2

如果你有指数衰减，基础离散概率分布是 geometric distribution. (It's the discrete counterpart of the continuous exponential distribution。）这样的几何分布使用参数 p 与一次试验的成功概率（如有偏掷硬币）。该分布描述了获得一次成功所需的试验次数。

分布的预期均值是 1/p。因此，我们可以计算观察值的平均值来估计 p.

该函数作为 scipy.stats.geom 构成 scipy 的一部分。要对分布进行抽样，请使用 geom.rvs(estimated_p, size=2000).

下面是一些演示该方法的代码：

from scipy.stats import geom
import matplotlib.pyplot as plt
import numpy as np
from collections import defaultdict

observation_index = [1, 2, 3, 4, 7, 13]
observation_count = [352, 28, 8, 4, 4, 4]

observed_mean = sum([i * c for i, c in zip(observation_index, observation_count)]) / sum(observation_count)

estimated_p = 1 / observed_mean
print('observed_mean:', observed_mean)
print('estimated p:', estimated_p)

generated_values = geom.rvs(estimated_p, size=2000)
generated_dict = defaultdict(int)
for v in generated_values:
    generated_dict[v] += 1
generated_index = sorted(list (generated_dict.keys()))
generated_count = [generated_dict [i] for i in  generated_index]
print(generated_index)
print(generated_count)

输出：

observed_mean: 1.32
estimated p: 0.7575757575757576
new random sample:
    [1, 2, 3, 4, 5, 7]
    [1516, 365, 86, 26, 6, 1]

从概率分布生成随机变量

Generate random variables from a probability distribution

python

numpy

math

random

normal-distribution