基于频率数据有效地拟合 Scipy 中的分布答案

【问题标题】：Fitting Distributions in Scipy Based on Frequency Data Efficiently基于频率数据有效地拟合 Scipy 中的分布
【发布时间】：2021-04-05 20:53:35
【问题描述】：

我有一些数据想要适合分布。数据由频率给出。我的意思是，我有我观察到的每一个事件以及我观察到它的次数。所以像：

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 或其他类似工具中，是否可以根据频率拟合分布？如果没有，是否有更好的方法来拟合分布，而无需创建潜在的巨大值列表？

【问题讨论】：

寻找参数最常用的方法是最大似然法，在这种情况下，使用频率而不是单个数据完全等同于在每个数据上加上一个等于频率的权重.因此，您可以尝试寻找允许权重与拟合数据相关联的函数。我不知道 Scipy 是否允许这样做，也许它已经允许了。如果不是，也可以看看 R。做不到这一点，从头开始写也没什么大不了的。
谢谢@RobertDodier。似乎 Scipy 不允许使用权重。

标签： python scipy statistics

【解决方案1】：

不幸的是，看着source，数据的“物化”方面似乎是硬编码的。不过，该功能并不复杂，因此您可以制作自己的版本。 TBH，如果你的总 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")

给我

【讨论】：

感谢您提供的信息。我不知道石榴。对于我的问题，扩展到样本会创建一个包含大约 10,000 个值的数组。我也在对几十万个样本进行拟合。使用 Pomegranate 将我的运行时间从大约 35 分钟缩短到大约 2 分钟。很酷！

【解决方案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)

【讨论】：