我有一个对数正态的分布式样本集 . 我可以使用具有线性或对数x轴的组织图来可视化样本 . 我可以对直方图进行拟合以获得PDF,然后将其缩放到具有线性x轴的图中的histrogram,另请参见this previously posted question .
但是,我无法使用对数x轴将PDF正确绘制到绘图中 .
不幸的是,将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()
Update 1 :
我忘了提到我使用的版本:
python 2.7.6
numpy 1.8.2
matplotlib 1.3.1
scipy 0.13.3
Update 2 :
正如@Christoph和@zaxliu所指出的那样(感谢两者),问题在于缩放PDF . 当我使用与直方图相同的箱子时,它可以工作,就像@ 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)
3 回答
根据我在@Warren Weckesser的原始答案中的理解,你reffered to "all you need to do"是:
我们可以尝试遵循他的建议,并根据箱子的中心点绘制两种获取pdf值的方法:
具有PDF功能
具有CDF功能:
您可以看到,第一个(使用pdf)和第二个(使用cdf)方法提供的结果几乎相同,并且两者都不完全匹配使用bin边缘计算的pdf .
如果放大,您会清楚地看到差异:
现在可以问的问题是:使用哪一个?我想答案将取决于但如果我们看一下累积概率:
您可以从输出中查看哪个方法更接近1.0:
CDF似乎给出了最接近的近似值 .
这很长,但我希望这是有道理的 .
Update:
我已经调整了代码来说明如何平滑PDF行 . 注意
s
变量,它定义了线的平滑程度 . 我在变量中添加了_s
后缀,以指示调整需要发生的位置 .这产生了这个情节:
如果放大平滑版本与非平滑版本,您会看到:
希望这可以帮助 .
由于我遇到了同样的问题并弄明白了,我想解释一下是什么,并为原始问题提供不同的解决方案 .
当您使用对数箱进行直方图时,这相当于更改变量
,其中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并除以上面计算的面积:
导致以下情节:
或者,您可以通过在日志空间中集成直方图来更改直方图“区域”的定义 . 请记住,在日志空间(t变量)中,PDF具有区域1.因此,您可以跳过缩放因子,并将上面的行替换为:
后一种解决方案可能是优选的,因为它不依赖于有关手头分布的任何信息 . 它适用于任何分发,而不仅仅是log-normal .
作为参考,这是添加了我的行的原始脚本:
正如@Christoph指出的那样,问题在于你缩放采样pdf的方式 .
因为pdf是概率密度的密度,如果你想要一个bin中的预期频率,你应该首先将密度乘以bin长度得到近似值样本将落入此区间的概率,然后您可以将此概率乘以样本总数,以估计将落入此区域的样本数 .
换句话说,每个bin应该以对数标度不均匀地缩放,而你用“hist下的区域”统一缩放它们 . 作为修复,您可以执行以下操作:
此外,您可能还需要考虑以类似的方式修改线性比例的缩放方法 . 实际上,您不需要累积面积,只需要按容器大小和样本总数计算多个密度 .
更新
在我看来,我目前估算箱子中概率的方法可能不是最准确的 . 由于pdf曲线是凹的,因此使用中点上的样本进行估计可能更准确 .