使用scipy拟合给定直方图的分布

2024-03-29 14:54:10 发布

您现在位置:Python中文网/ 问答频道 /正文

我想用scipy(在我的例子中,使用weibull_min)来拟合我的数据分布。在给定直方图而不是数据点的情况下,有没有可能做到这一点?在我的例子中,因为直方图有大小为1的整数存储箱,我知道我可以用以下方式推断数据:

import numpy as np
orig_hist = np.array([10, 5, 3, 2, 1])

ext_data = reduce(lambda x,y: x+y, [[i]*x for i, x in enumerate(orig_hist)])

在这种情况下,ext_data将保存:

^{pr2}$

并使用以下方法构建直方图:

np.histogram(ext_data, bins=5)

相当于原始历史

然而,鉴于我已经建立了直方图,我希望避免外推数据,并使用orig峎hist来拟合分布,但我不知道是否可以在拟合过程中直接使用它。另外,是否有一个numpy函数可以用来执行类似于我所展示的推断的功能?在


Tags: 数据numpydatanp情况整数scipy直方图
1条回答
网友
1楼 · 发布于 2024-03-29 14:54:10

我可能误解了一些事情,但我相信,拟合直方图正是您应该做的:您试图近似概率密度。最接近的概率是直方图。你只需要规范化它就可以得到1的积分,或者允许你的模型包含一个任意的前置因子。在

import numpy as np
import scipy.stats as stats
import scipy.optimize as opt
import matplotlib.pyplot as plt

orig_hist = np.array([10, 5, 3, 2, 1])
norm_hist = orig_hist/float(sum(orig_hist))

popt,pcov = opt.curve_fit(lambda x,c: stats.weibull_min.pdf(x,c), np.arange(len(norm_hist)),norm_hist)

plt.figure()
plt.plot(norm_hist,'o-',label='norm_hist')
plt.plot(stats.weibull_min.pdf(np.arange(len(norm_hist)),popt),'s-',label='Weibull_min fit')
plt.legend()

当然,对于给定的输入,Weibull拟合远远不能令人满意:

fit to data

更新

如前所述,Weibull-min不适合您的样本输入。更大的问题是,它与实际数据的拟合度也很差:

^{pr2}$

new histogram data

这个直方图有两个主要问题。第一,正如我所说,它不太可能对应于威布尔分布:它是接近零的最大分布,并且有一个长尾,因此它需要威布尔参数的非平凡组合。此外,你的直方图显然只包含了分布的一部分。这意味着我上面的规范化建议肯定会失败。您无法避免在拟合中使用任意比例参数。在

我手动定义了一个标度Weibull拟合函数according to the formula on Wikipedia

my_weibull = lambda x,l,c,A: A*float(c)/l*(x/float(l))**(c-1)*np.exp(-(x/float(l))**c)

在这个函数中,x是自变量,llambda(缩放参数),c是{}(形状参数),A是缩放前置因子。引入A的一个微弱的好处是,您不必对直方图进行标准化。在

现在,当我把这个函数放到scipy.optimize.curve_fit中时,我发现了您所做的:它实际上并不执行拟合,而是坚持初始拟合参数,无论您设置什么(使用p0参数;每个参数的默认猜测都是1)。curve_fit似乎认为拟合收敛。在

在一个多小时的与墙壁有关的头部撞击之后,我意识到问题在于x=0处的奇异行为导致了非线性最小二乘算法的失效。通过排除您的第一个数据点,您就获得了与数据的实际拟合。我怀疑,如果我们设置c=1并且不允许它适合,那么这个问题可能会消失,但是允许它被安装可能更有用(所以我没有检查)。在

下面是相应的代码:

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

orig_hist = np.array([ 23., 14., 13., 12., 12., 12., 11., 11., 11., 11., 10., 10., 10., 9., 9., 8., 8., 8., 8., 8., 8., 8., 8., 8., 8., 8., 7., 7., 7., 7., 7., 7., 7., 7., 7., 7., 7., 7., 7., 6., 6., 6., 6., 6., 6., 6., 6., 6., 6., 6.], dtype=np.float32)

my_weibull = lambda x,l,c,A: A*float(c)/l*(x/float(l))**(c-1)*np.exp(-(x/float(l))**c)

popt,pcov = opt.curve_fit(my_weibull,np.arange(len(orig_hist))[1:],orig_hist[1:]) #throw away x=0!

plt.figure()
plt.plot(np.arange(len(orig_hist)),orig_hist,'o-',label='orig_hist')
plt.plot(np.arange(len(orig_hist)),my_weibull(np.arange(len(orig_hist)),*popt),'s-',label='Scaled Weibull fit')
plt.legend()

结果:

new fit

In [631]: popt
Out[631]: array([  1.10511850e+02,   8.82327822e-01,   1.05206207e+03])

最终拟合的参数顺序为(l,c,A),形状参数约为0.88。这对应于一个发散的概率密度,这解释了为什么会出现一些错误

RuntimeWarning: invalid value encountered in power

以及为什么没有来自x=0拟合的数据点。但是从数据和拟合度之间的视觉一致性来看,您可以评估结果是否可以接受。在

如果你想做得过火,你可以尝试用这些参数使用np.random.weibull生成点,然后将得到的直方图与你自己的直方图进行比较。在

相关问题 更多 >