在 python 中使用对数轴缩放和拟合对数正态分布

Question

我有一组对数正态分布的样本。我可以使用具有线性或对数 x 轴的直方图来可视化样本。我可以对直方图进行拟合以获取 PDF，然后使用线性 x 轴将其缩放到图中的直方图，另请参见 this previously posted question。

但是，我无法正确地将 PDF 绘制到带有对数 x 轴的图中。

不幸的是，这不仅是 PDF 区域缩放到直方图的问题，而且 PDF 也会向左移动，如下图所示。

我现在的问题是，我在这里做错了什么？使用 CDF 绘制预期的直方图 as suggested in this answer 是可行的。我只是想知道我在这段代码中做错了什么，因为据我所知它也应该有效。

这是python代码（对不起，它比较长，但我想发布一个“完整的单机版”）：

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

# generate log-normal distributed set of samples
np.random.seed(42)
samples   = np.random.lognormal( mean=1, sigma=.4, size=10000 )

# make a fit to the samples
shape, loc, scale = scipy.stats.lognorm.fit( samples, floc=0 )
x_fit       = np.linspace( samples.min(), samples.max(), 100 )
samples_fit = scipy.stats.lognorm.pdf( x_fit, shape, loc=loc, scale=scale )

# plot a histrogram with linear x-axis
plt.subplot( 1, 2, 1 )
N_bins = 50
counts, bin_edges, ignored = plt.hist( samples, N_bins, histtype='stepfilled', label='histogram' )
# calculate area of histogram (area under PDF should be 1)
area_hist = .0
for ii in range( counts.size):
    area_hist += (bin_edges[ii+1]-bin_edges[ii]) * counts[ii]
# oplot fit into histogram
plt.plot( x_fit, samples_fit*area_hist, label='fitted and area-scaled PDF', linewidth=2)
plt.legend()

# make a histrogram with a log10 x-axis
plt.subplot( 1, 2, 2 )
# equally sized bins (in log10-scale)
bins_log10 = np.logspace( np.log10( samples.min()  ), np.log10( samples.max() ), N_bins )
counts, bin_edges, ignored = plt.hist( samples, bins_log10, histtype='stepfilled', label='histogram' )
# calculate area of histogram
area_hist_log = .0
for ii in range( counts.size):
    area_hist_log += (bin_edges[ii+1]-bin_edges[ii]) * counts[ii]
# get pdf-values for log10 - spaced intervals
x_fit_log       = np.logspace( np.log10( samples.min()), np.log10( samples.max()), 100 )
samples_fit_log = scipy.stats.lognorm.pdf( x_fit_log, shape, loc=loc, scale=scale )
# oplot fit into histogram
plt.plot( x_fit_log, samples_fit_log*area_hist_log, label='fitted and area-scaled PDF', linewidth=2 )

plt.xscale( 'log' )
plt.xlim( bin_edges.min(), bin_edges.max() )
plt.legend()
plt.show()

更新 1：

我忘了提及我正在使用的版本：

python      2.7.6
numpy       1.8.2
matplotlib  1.3.1
scipy       0.13.3

更新 2：

正如@Christoph 和@zaxliu 所指出的（感谢两者），问题在于PDF 的缩放。当我使用与直方图相同的 bin 时，它可以工作，就像在@zaxliu 的解决方案中一样，但是在为 PDF 使用更高分辨率时我仍然遇到一些问题（如我上面的示例所示）。如下图所示：

右侧图的代码是（我省略了导入和数据样本生成的东西，你可以在上面的例子中找到它们）：

# equally sized bins in log10-scale
bins_log10 = np.logspace( np.log10( samples.min()  ), np.log10( samples.max() ), N_bins )
counts, bin_edges, ignored = plt.hist( samples, bins_log10, histtype='stepfilled', label='histogram' )

# calculate length of each bin (required for scaling PDF to histogram)
bins_log_len = np.zeros( bins_log10.size )
for ii in range( counts.size):
    bins_log_len[ii] = bin_edges[ii+1]-bin_edges[ii]

# get pdf-values for same intervals as histogram
samples_fit_log = scipy.stats.lognorm.pdf( bins_log10, shape, loc=loc, scale=scale )

# oplot fitted and scaled PDF into histogram
plt.plot( bins_log10, np.multiply(samples_fit_log,bins_log_len)*sum(counts), label='PDF using histogram bins', linewidth=2 )

# make another pdf with a finer resolution
x_fit_log       = np.logspace( np.log10( samples.min()), np.log10( samples.max()), 100 )
samples_fit_log = scipy.stats.lognorm.pdf( x_fit_log, shape, loc=loc, scale=scale )
# calculate length of each bin (required for scaling PDF to histogram)
# in addition, estimate middle point for more accuracy (should in principle also be done for the other PDF)
bins_log_len       = np.diff( x_fit_log )
samples_log_center = np.zeros( x_fit_log.size-1 )
for ii in range( x_fit_log.size-1 ):
    samples_log_center[ii] = .5*(samples_fit_log[ii] + samples_fit_log[ii+1] )

# scale PDF to histogram
# NOTE: THIS IS NOT WORKING PROPERLY (SEE FIGURE)
pdf_scaled2hist = np.multiply(samples_log_center,bins_log_len)*sum(counts)

# oplot fit into histogram
plt.plot( .5*(x_fit_log[:-1]+x_fit_log[1:]), pdf_scaled2hist, label='PDF using own bins', linewidth=2 )

plt.xscale( 'log' )
plt.xlim( bin_edges.min(), bin_edges.max() )
plt.legend(loc=3)

Answer 1

正如@Christoph 所指出的，问题在于您缩放采样 pdf 的方式。

因为 pdf 是概率密度的密度，如果你想要一个 bin 中的预期频率，你应该先将密度乘以 bin 长度，得到一个样本落入这个 bin 的近似概率，然后你可以乘以这个概率通过样本总数来估计将落入这个 bin 的样本数。

换句话说，每个 bin 应该以对数尺度不均匀地缩放，而您使用“hist under hist”统一缩放它们。作为修复，您可以执行以下操作：

# make a histrogram with a log10 x-axis
plt.subplot( 1, 2, 2 )
# equally sized bins (in log10-scale)
bins_log10 = np.logspace( np.log10( samples.min()  ), np.log10( samples.max() ), N_bins )
counts, bin_edges, ignored = plt.hist( samples, bins_log10, histtype='stepfilled', label='histogram' )
# calculate length of each bin
len_bin_log = np.zeros([bins_log10.size,])
for ii in range( counts.size):
    len_bin_log[ii] = (bin_edges[ii+1]-bin_edges[ii])

# get pdf-values for log10 - spaced intervals
# x_fit_log       = np.logspace( np.log10( samples.min()), np.log10( samples.max()), N_bins )
samples_fit_log = scipy.stats.lognorm.pdf( bins_log10, shape, loc=loc, scale=scale )

# oplot fit into histogram
plt.plot(bins_log10 , np.multiply(samples_fit_log,len_bin_log)*sum(counts), label='fitted and area-scaled PDF', linewidth=2 )
plt.xscale( 'log' )
plt.xlim( bin_edges.min(), bin_edges.max() )
# plt.legend()
plt.show()

此外，您可能还想考虑以类似的方式修改线性比例的缩放方法。实际上，您不需要累积面积，只需按 bin 大小乘以密度乘以样本总数。

更新

我想到，我目前估计垃圾箱概率的方法可能不是最准确的方法。由于 pdf 曲线是凹的，因此在中间点使用样本进行估计可能会更准确。

Answer 2

根据我在@Warren Weckesser 的原始答案中的理解，您reffered to“所有你需要做的”是：

将cdf(b) - cdf(a) 的近似值写为 cdf(b) - cdf(a) = pdf(m)*(b - a) 其中，m 是区间 [a, b] 的中点

我们可以尝试按照他的建议，根据 bin 的中心点绘制两种获取 pdf 值的方法：

1. 带PDF功能
2. 带 CDF 功能：

---

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats 

# generate log-normal distributed set of samples
np.random.seed(42)
samples = np.random.lognormal(mean=1, sigma=.4, size=10000)
N_bins = 50

# make a fit to the samples
shape, loc, scale = stats.lognorm.fit(samples, floc=0)
x_fit       = np.linspace(samples.min(), samples.max(), 100)
samples_fit = stats.lognorm.pdf(x_fit, shape, loc=loc, scale=scale)

# plot a histrogram with linear x-axis
fig, (ax1, ax2) = plt.subplots(1,2, figsize=(10,5), gridspec_kw={'wspace':0.2})
counts, bin_edges, ignored = ax1.hist(samples, N_bins, histtype='stepfilled', alpha=0.4,
                                      label='histogram')

# calculate area of histogram (area under PDF should be 1)
area_hist = ((bin_edges[1:] - bin_edges[:-1]) * counts).sum()

# plot fit into histogram
ax1.plot(x_fit, samples_fit*area_hist, label='fitted and area-scaled PDF', linewidth=2)
ax1.legend()

# equally sized bins in log10-scale and centers
bins_log10 = np.logspace(np.log10(samples.min()), np.log10(samples.max()), N_bins)
bins_log10_cntr = (bins_log10[1:] + bins_log10[:-1]) / 2

# histogram plot
counts, bin_edges, ignored = ax2.hist(samples, bins_log10, histtype='stepfilled', alpha=0.4,
                                      label='histogram')

# calculate length of each bin and its centers(required for scaling PDF to histogram)
bins_log_len = np.r_[bin_edges[1:] - bin_edges[: -1], 0]
bins_log_cntr = bin_edges[1:] - bin_edges[:-1]

# get pdf-values for same intervals as histogram
samples_fit_log = stats.lognorm.pdf(bins_log10, shape, loc=loc, scale=scale)

# pdf-values for centered scale
samples_fit_log_cntr = stats.lognorm.pdf(bins_log10_cntr, shape, loc=loc, scale=scale)

# pdf-values using cdf 
samples_fit_log_cntr2_ = stats.lognorm.cdf(bins_log10, shape, loc=loc, scale=scale)
samples_fit_log_cntr2 = np.diff(samples_fit_log_cntr2_)

# plot fitted and scaled PDFs into histogram
ax2.plot(bins_log10, 
         samples_fit_log * bins_log_len * counts.sum(), '-', 
         label='PDF with edges',  linewidth=2)

ax2.plot(bins_log10_cntr, 
         samples_fit_log_cntr * bins_log_cntr * counts.sum(), '-', 
         label='PDF with centers', linewidth=2)

ax2.plot(bins_log10_cntr, 
         samples_fit_log_cntr2 * counts.sum(), 'b-.', 
         label='CDF with centers', linewidth=2)


ax2.set_xscale('log')
ax2.set_xlim(bin_edges.min(), bin_edges.max())
ax2.legend(loc=3)
plt.show()

您可以看到，第一种（使用 pdf）和第二种（使用 cdf）方法都给出了几乎相同的结果，并且两者都不完全匹配使用 bin 边缘计算的 pdf。

如果放大，您会清楚地看到差异：

现在人们可以提出的问题是：使用哪一个？我想答案将取决于但如果我们看一下累积概率：

print 'Cumulative probabilities:'
print 'Using edges:         {:>10.5f}'.format((samples_fit_log * bins_log_len).sum())
print 'Using PDF of centers:{:>10.5f}'.format((samples_fit_log_cntr * bins_log_cntr).sum())
print 'Using CDF of centers:{:>10.5f}'.format(samples_fit_log_cntr2.sum())

您可以从输出中看到哪个方法更接近 1.0：

Cumulative probabilities:
Using edges:            1.03263
Using PDF of centers:   0.99957
Using CDF of centers:   0.99991

CDF 似乎给出了最接近的近似值。

这很长，但我希望这是有道理的。

更新：

我已经调整了代码来说明如何平滑 PDF 线条。 注意s 变量，它定义了线条的平滑程度。 我在变量中添加了_s 后缀，以指示需要在哪里进行调整。

# generate log-normal distributed set of samples
np.random.seed(42)
samples = np.random.lognormal(mean=1, sigma=.4, size=10000)
N_bins = 50

# make a fit to the samples
shape, loc, scale = stats.lognorm.fit(samples, floc=0)

# plot a histrogram with linear x-axis
fig, ax2 = plt.subplots()#1,2, figsize=(10,5), gridspec_kw={'wspace':0.2})

# equally sized bins in log10-scale and centers
bins_log10 = np.logspace(np.log10(samples.min()), np.log10(samples.max()), N_bins)
bins_log10_cntr = (bins_log10[1:] + bins_log10[:-1]) / 2

# smoother PDF line
s = 10 # mulpiplier to N_bins - the bigger s is the smoother the line
bins_log10_s = np.logspace(np.log10(samples.min()), np.log10(samples.max()), N_bins * s)
bins_log10_cntr_s = (bins_log10_s[1:] + bins_log10_s[:-1]) / 2

# histogram plot
counts, bin_edges, ignored = ax2.hist(samples, bins_log10, histtype='stepfilled', alpha=0.4,
                                      label='histogram')

# calculate length of each bin and its centers(required for scaling PDF to histogram)
bins_log_len = np.r_[bins_log10_s[1:] - bins_log10_s[: -1], 0]
bins_log_cntr = bins_log10_s[1:] - bins_log10_s[:-1]

# smooth pdf-values for same intervals as histogram
samples_fit_log_s = stats.lognorm.pdf(bins_log10_s, shape, loc=loc, scale=scale)

# pdf-values for centered scale
samples_fit_log_cntr = stats.lognorm.pdf(bins_log10_cntr_s, shape, loc=loc, scale=scale)

# smooth pdf-values using cdf 
samples_fit_log_cntr2_s_ = stats.lognorm.cdf(bins_log10_s, shape, loc=loc, scale=scale)
samples_fit_log_cntr2_s = np.diff(samples_fit_log_cntr2_s_)

# plot fitted and scaled PDFs into histogram
ax2.plot(bins_log10_cntr_s, 
         samples_fit_log_cntr * bins_log_cntr * counts.sum() * s, '-', 
         label='Smooth PDF with centers', linewidth=2)

ax2.plot(bins_log10_cntr_s, 
         samples_fit_log_cntr2_s * counts.sum() * s, 'k-.', 
         label='Smooth CDF with centers', linewidth=2)

ax2.set_xscale('log')
ax2.set_xlim(bin_edges.min(), bin_edges.max())
ax2.legend(loc=3)
plt.show)

这会产生这个情节：

如果您放大平滑版本与非平滑版本，您会看到：

希望这会有所帮助。

Answer 3

由于我遇到了同样的问题并想通了，我想解释一下发生了什么，并为原始问题提供不同的解决方案。

当您使用对数 bin 绘制直方图时，这相当于更改变量 ，其中 x 是您的原始样本（或您用来绘制它们的网格），t 是相对于箱是线性间隔的。因此，实际对应直方图的PDF是

我们仍在使用 x 变量作为 PDF 的输入，所以这变成了

您需要将 PDF 乘以 x！

这固定了 PDF 的形状，但我们仍然需要缩放 PDF 以使曲线下的面积等于直方图。事实上，PDF 下的面积不等于 1，因为我们是在 x 上积分，并且

因为我们正在处理对数正态变量。因为根据scipy documentation，分布参数对应shape = sigma和scale = exp(mu)，我们可以很容易地计算出你代码中的右边为scale * np.exp(shape**2/2.)。

事实上，一行代码修复了您的原始脚本，将计算出的 PDF 值乘以 x 并除以上面计算的面积：

samples_fit_log *= x_fit_log / (scale * np.exp(shape**2/2.))

导致以下情节：

或者，您可以通过在日志空间中集成直方图来更改直方图“区域”的定义。请记住，在对数空间（t 变量）中，PDF 的面积为 1。因此您可以跳过缩放因子，并将上面的行替换为：

area_hist_log = np.dot(np.diff(np.log(bin_edges)), counts)
samples_fit_log *= x_fit_log

后一种解决方案可能更可取，因为它不依赖于任何有关手头分布的信息。它适用于任何分布，而不仅仅是对数正态分布。

作为参考，这是添加了我的行的原始脚本：

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

# generate log-normal distributed set of samples
np.random.seed(42)
samples   = np.random.lognormal( mean=1, sigma=.4, size=10000 )

# make a fit to the samples
shape, loc, scale = scipy.stats.lognorm.fit( samples, floc=0 )
x_fit       = np.linspace( samples.min(), samples.max(), 100 )
samples_fit = scipy.stats.lognorm.pdf( x_fit, shape, loc=loc, scale=scale )

# plot a histrogram with linear x-axis
plt.subplot( 1, 2, 1 )
N_bins = 50
counts, bin_edges, ignored = plt.hist( samples, N_bins, histtype='stepfilled', label='histogram' )
# calculate area of histogram (area under PDF should be 1)
area_hist = .0
for ii in range( counts.size):
    area_hist += (bin_edges[ii+1]-bin_edges[ii]) * counts[ii]
# oplot fit into histogram
plt.plot( x_fit, samples_fit*area_hist, label='fitted and area-scaled PDF', linewidth=2)
plt.legend()

# make a histrogram with a log10 x-axis
plt.subplot( 1, 2, 2 )
# equally sized bins (in log10-scale)
bins_log10 = np.logspace( np.log10( samples.min()  ), np.log10( samples.max() ), N_bins )
counts, bin_edges, ignored = plt.hist( samples, bins_log10, histtype='stepfilled', label='histogram' )
# calculate area of histogram
area_hist_log = .0
for ii in range( counts.size):
    area_hist_log += (bin_edges[ii+1]-bin_edges[ii]) * counts[ii]
# get pdf-values for log10 - spaced intervals
x_fit_log       = np.logspace( np.log10( samples.min()), np.log10( samples.max()), 100 )
samples_fit_log = scipy.stats.lognorm.pdf( x_fit_log, shape, loc=loc, scale=scale )
# scale pdf output:
samples_fit_log *= x_fit_log / (scale * np.exp(shape**2/2.))
# alternatively you could do:
#area_hist_log = np.dot(np.diff(np.log(bin_edges)), counts)
#samples_fit_log *= x_fit_log

# oplot fit into histogram
plt.plot( x_fit_log, samples_fit_log*area_hist_log, label='fitted and area-scaled PDF', linewidth=2 )

plt.xscale( 'log' )
plt.xlim( bin_edges.min(), bin_edges.max() )
plt.legend()
plt.show()