使用 scipy 拟合给定直方图的分布
Fitting a distribution given the histogram using scipy
我想使用 scipy（在我的例子中，使用 weibull_min）对我的数据进行分布拟合。给定直方图而不是数据点是否可以做到这一点？在我的例子中，因为直方图有大小为 1 的整数 bin，我知道我可以通过以下方式推断我的数据：
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 将持有：
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4]
并使用以下方法构建直方图：
np.histogram(ext_data, bins=5)
等同于 orig_hist
然而，鉴于我已经构建了直方图，我想避免外推数据并使用 orig_hist 来拟合分布，但我不知道是否可以直接使用它在装修过程中。此外，是否有一个 numpy 函数可用于执行与我展示的外推类似的操作？
我可能误解了什么，但我相信拟合直方图正是您应该做的：您正在尝试近似概率密度。直方图尽可能接近潜在的概率密度。您只需对其进行归一化以获得 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 拟合将远不能令人满意：
更新
正如我上面提到的，Weibull_min 不适合您的样本输入。更大的问题是它也不太适合你的实际数据：
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)
这个直方图有两个主要问题。首先，正如我所说，它不太可能对应于 Weibull_min 分布：它在零附近达到最大值并且尾巴很长，因此它需要 Weibull 参数的 non-trivial 组合。此外，您的直方图显然只包含分布的一部分。这意味着我上面的规范化建议肯定会失败。您无法避免在拟合中使用任意比例参数。
我手动定义了一个缩放的威布尔拟合函数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是自变量，l是lambda（尺度参数），c是k（形状参数) 并且 A 是缩放前因子。引入 A 的微弱好处是您不必标准化直方图。
现在，当我将这个函数放到 scipy.optimize.curve_fit 中时，我发现了你所做的：它实际上并没有执行拟合，而是坚持使用初始拟合参数，无论你设置什么（使用 p0 参数；每个参数的默认猜测都是 1）。 和 curve_fit好像认为拟合收敛了
经过一个多小时的 wall-related head-banging，我意识到问题是 x=0 处的奇异行为抛出了非线性 least-squares 算法。通过排除您的 very first 数据点，您可以实际拟合您的数据。我怀疑如果我们设置 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()
结果：
In [631]: popt
Out[631]: array([ 1.10511850e+02, 8.82327822e-01, 1.05206207e+03])
最终拟合的参数顺序(l,c,A)，shape参数在0.88左右。这对应于发散的概率密度，这解释了为什么弹出一些错误说
RuntimeWarning: invalid value encountered in power
以及为什么 x=0 的拟合没有数据点。但是从数据和拟合的视觉一致性来看，你可以评估结果是否可以接受。
如果你想过度，你可以尝试使用 np.random.weibull 和这些参数生成点，然后将生成的直方图与你自己的直方图进行比较。
我想使用 scipy（在我的例子中，使用 weibull_min）对我的数据进行分布拟合。给定直方图而不是数据点是否可以做到这一点？在我的例子中，因为直方图有大小为 1 的整数 bin，我知道我可以通过以下方式推断我的数据：
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 将持有：
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4]
并使用以下方法构建直方图：
np.histogram(ext_data, bins=5)
等同于 orig_hist
然而，鉴于我已经构建了直方图，我想避免外推数据并使用 orig_hist 来拟合分布，但我不知道是否可以直接使用它在装修过程中。此外，是否有一个 numpy 函数可用于执行与我展示的外推类似的操作？
我可能误解了什么，但我相信拟合直方图正是您应该做的：您正在尝试近似概率密度。直方图尽可能接近潜在的概率密度。您只需对其进行归一化以获得 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 拟合将远不能令人满意：
更新
正如我上面提到的，Weibull_min 不适合您的样本输入。更大的问题是它也不太适合你的实际数据：
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)
这个直方图有两个主要问题。首先，正如我所说，它不太可能对应于 Weibull_min 分布：它在零附近达到最大值并且尾巴很长，因此它需要 Weibull 参数的 non-trivial 组合。此外，您的直方图显然只包含分布的一部分。这意味着我上面的规范化建议肯定会失败。您无法避免在拟合中使用任意比例参数。
我手动定义了一个缩放的威布尔拟合函数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是自变量，l是lambda（尺度参数），c是k（形状参数) 并且 A 是缩放前因子。引入 A 的微弱好处是您不必标准化直方图。
现在，当我将这个函数放到 scipy.optimize.curve_fit 中时，我发现了你所做的：它实际上并没有执行拟合，而是坚持使用初始拟合参数，无论你设置什么（使用 p0 参数；每个参数的默认猜测都是 1）。 和 curve_fit好像认为拟合收敛了
经过一个多小时的 wall-related head-banging，我意识到问题是 x=0 处的奇异行为抛出了非线性 least-squares 算法。通过排除您的 very first 数据点，您可以实际拟合您的数据。我怀疑如果我们设置 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()
结果：
In [631]: popt
Out[631]: array([ 1.10511850e+02, 8.82327822e-01, 1.05206207e+03])
最终拟合的参数顺序(l,c,A)，shape参数在0.88左右。这对应于发散的概率密度，这解释了为什么弹出一些错误说
RuntimeWarning: invalid value encountered in power
以及为什么 x=0 的拟合没有数据点。但是从数据和拟合的视觉一致性来看，你可以评估结果是否可以接受。
如果你想过度，你可以尝试使用 np.random.weibull 和这些参数生成点，然后将生成的直方图与你自己的直方图进行比较。