在Scipy中基于频率数据高效地拟合分布

Question

我有一些数据需要拟合到一个分布中。这些数据是按频率给出的。我的意思是，我记录了每个事件和我观察到它的次数。例如：

data = [(1, 34), (2, 1023), (3, 3243), (4, 879), (5, 202), (6, 10)]

每个元组中的第一个数字是我观察到的事件，第二个数字是该事件的总观察次数。

使用Scipy，我可以通过调用scipy.stats.lognorm.fit来拟合（例如）对数正态分布。然而，这个例程期望看到所有观察结果的列表，而不是频率。我可以像这样拟合分布：

import scipy
temp_data = []
for x in data:
    temp_data += [x[0]] * x[1]
params = scipy.stats.lognorm.fit(temp_data)

但是，哇，那看起来非常低效。

在Scipy或其他类似工具中，是否有一种基于频率拟合分布的方法？如果没有，是否有更好的方法可以拟合分布而不必创建一个可能巨大的值列表？

- Neither_8

---

找到参数的最常见方法是最大似然法，这种情况下，使用频率而不是个体数据进行计算，等同于给每个数据赋予权重，该权重等于频率。因此，您可以尝试寻找允许在拟合过程中将权重与数据相关联的函数。我不知道Scipy是否允许这样做，也许它已经支持了。如果不行，您可以考虑使用R语言。如果都不行，从头开始编写也不是什么大问题。 - Robert Dodier

感谢 @RobertDodier。看起来 Scipy 不允许使用权重。 - Neither_8

Answer 1

很不幸，看起来来源中的数据“物化”方面是硬编码的。虽然该函数并不那么复杂，但你可以自己制作一个版本。老实说，如果你的总N仍然可控，我可能会只做data = np.array(data); expanded_data = np.repeat(data[:,0], data[:,1])尽管效率低下，因为生命苦短。

另一种选择是使用pomegranate，它支持传递权重：

import numpy as np
import scipy.stats
import matplotlib.pyplot as plt
import pomegranate as pg

data = [(1, 34), (2, 1023), (3, 3243), (4, 879), (5, 202), (6, 10)]

data = np.array(data)
expanded = np.repeat(data[:,0], data[:,1].astype(int))

scipy_shape, _, scipy_scale = scipy_params = scipy.stats.lognorm.fit(expanded, floc=0)
scipy_sigma, scipy_mu = scipy_shape, np.log(scipy_scale)

pg_dist = pg.LogNormalDistribution(0, 1)
pg_dist.fit(data[:,0], weights=data[:,1])
pg_mu, pg_sigma = pg_dist.parameters

fig = plt.figure()
ax = fig.add_subplot(111)

x = np.linspace(0.1, 10, 100)
ax.plot(data[:,0], data[:, 1] / data[:,1].sum(), label="freq")
ax.plot(x, scipy.stats.lognorm(*scipy_params).pdf(x),
        label=r"scipy: $\mu$ {:1.3f} $\sigma$ {:1.3f}".format(scipy_mu, scipy_sigma), alpha=0.5)
ax.plot(x, pg_dist.probability(x),
        label=r"pomegranate: $\mu$ {:1.3f} $\sigma$ {:1.3f}".format(pg_mu, pg_sigma), linestyle='--', alpha=0.5)
ax.legend(loc='upper right')
fig.savefig("compare.png")

给我

- DSM

---

谢谢你提供的信息。我之前不知道有pomegranate这个工具。对于我的问题，将样本扩展成数组后大约有10,000个值。而且我还要对几十万个样本集进行拟合。使用Pomegranate可以将运行时间从大约35分钟缩短到2分钟左右。非常棒！ - Neither_8

Answer 2

你可以根据频率分布随机抽取样本，并进行拟合：

import scipy
import numpy as np

data = np.array(
    [(1, 34), (2, 1023), (3, 3243), (4, 879), (5, 202), (6, 10)], 
    dtype=float,
)
values = data[0]
weights = data[1]
seed = 87

gen = np.random.default_rng(seed)
sample = gen.choices(
    values, size=500, p=weights/np.sum(weights))

params = scipy.stats.lognorm.fit(values)

- shaunc