将对数正态分布拟合到已分箱的数据

Question

我想对我已经分箱的数据进行对数正态拟合。条形图如下所示：

不幸的是，当我尝试使用标准

lognorm.pdf()

时，拟合分布的形状非常不同。我想这是因为我的数据已经被装箱了。这是代码：

times, data, bin_points = ReadHistogramFile(filename)

xmin = 200
xmax = 800
x = np.linspace(xmin, xmax, 1000)
shape, loc, scale = stats.lognorm.fit(data, floc=0)
pdf = stats.lognorm.pdf(x, shape, loc=loc, scale=scale)

area=data.sum()
plt.bar(bars, data, width=10, color='b')
plt.plot(x*area, pdf, 'k' )

拟合分布如下： 显然缩放也有问题。不过我对此不太关心。我的主要问题是分布的形状。这可能与：这个问题重复，但我找不到正确的解决方案。我尝试了一下，仍然得到了与执行上述操作时非常相似的形状。感谢您的帮助！

更新： 通过使用

curve_fit()

，我能够稍微适应一下。但我还不满意。我想要原来的垃圾箱而不是统一垃圾箱。我也不确定到底发生了什么，以及是否没有更好的选择。这是代码：

def normalize_integral(data, bin_size):
normalized_data = np.zeros(size(data))
print bin_size
sum = data.sum()
integral = bin_size*sum
for i in range(0, size(data)-1):
    normalized_data[i] = data[i]/integral
    
print 'integral:', normalized_data.sum()*bin_size
return normalized_data



def pdf(x, mu, sigma):
"""pdf of lognormal distribution"""

return (np.exp(-(np.log(x) - mu)**2 / (2 * sigma**2)) / (x * sigma * np.sqrt(2 * np.pi)))


bin_points=np.linspace(280.5, 1099.55994, len(bin_points))
data=[9.78200000e+03 1.15120000e+04 1.18000000e+04 1.79620000e+04 2.76980000e+04   2.78260000e+04   3.35460000e+04   3.24260000e+04 3.16500000e+04   3.30820000e+04   4.84560000e+04   5.86500000e+04 6.34220000e+04   5.11880000e+04   5.13180000e+04   4.74320000e+04 4.35420000e+04   4.13400000e+04   3.60880000e+04   2.96900000e+04 2.66640000e+04   2.58720000e+04   2.57560000e+04   2.20960000e+04 1.46880000e+04   9.97200000e+03   5.74200000e+03   3.52000000e+03 2.74600000e+03   2.61800000e+03   1.50000000e+03   7.96000000e+02 5.40000000e+02   2.98000000e+02   2.90000000e+02   2.22000000e+02 2.26000000e+02   1.88000000e+02   1.20000000e+02   5.00000000e+01 5.40000000e+01   5.80000000e+01   5.20000000e+01   2.00000000e+01 2.80000000e+01   6.00000000e+00   0.00000000e+00   0.00000000e+00 0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00 0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00 0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00 0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00 0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00 0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00 0.00000000e+00   0.00000000e+00]
normalized_data_unitybins = normalize_integral(data,1)


plt.figure(figsize=(9,4))
ax1=plt.subplot(121)
ax2=plt.subplot(122)
ax2.bar(unity_bins, normalized_data_unitybins, width=1, color='b')
fitParams, fitCov = curve_fit(pdf, unity_bins, normalized_data_unitybins, p0=[1,1],maxfev = 1000000)
fitData=pdf(unity_bins, *fitParams)
ax2.plot(unity_bins, fitData,'g-')

ax1.bar(bin_points, normalized_data_unitybins, width=10, color='b')
fitParams, fitCov = curve_fit(pdf, bin_points, normalized_data_unitybins, p0=[1,1],maxfev = 1000000)
fitData=pdf(bin_points, *fitParams)
ax1.plot(bin_points, fitData,'g-')

Answer 1

正如您所提到的，您不能在分箱数据上使用

lognorm.fit

。所以你需要做的就是从直方图中恢复原始数据。显然这不是“无损”，垃圾箱越多越好。

带有一些生成数据的示例代码：

import numpy as np
import scipy.stats as stats
import matplotlib.pylab as plt

# generate some data
ln = stats.lognorm(0.4, scale=100)
data = ln.rvs(size=2000)

counts, bins, _ = plt.hist(data, bins=50, ec='k')
# note that the len of bins is 51, since it contains upper and lower limit of every bin

# restore data from histogram: counts multiplied bin centers
restored = [[d]*int(counts[n]) for n,d in enumerate((bins[1:]+bins[:-1])/2)]
# flatten the result
restored = [item for sublist in restored for item in sublist]

print(stats.lognorm.fit(restored, floc=0))

dist = stats.lognorm(*stats.lognorm.fit(restored, floc=0))
x = np.arange(1,400)
y = dist.pdf(x)

# the pdf is normalized, so we need to scale it to match the histogram
y = y/y.max()
y = y*counts.max()

plt.plot(x,y,'r',linewidth=2)
plt.show()