用对数拉普拉斯概率密度函数制作Matplotlib直方图

Question

(在深入阅读源代码之前，请务必查看本文末尾的编辑部分)

我正在绘制一个人口的直方图，它似乎是对数拉普拉斯分布：

我试图画出最佳拟合线来验证我的假设，但我无法得到有意义的结果。

我使用维基百科上的拉普拉斯概率密度函数定义，并取10的PDF幂（以“反转”对数直方图的效果）。

我做错了什么吗？

这是我的代码。我通过标准输入将数据传输 (cat pop.txt | python hist.py) -- 这里 有一个样本人口。

from pylab import *
import numpy    
def laplace(x, mu, b):
    return 10**(1.0/(2*b) * numpy.exp(-abs(x - mu)/b))    
def main():
    import sys
    num = map(int, sys.stdin.read().strip().split(' '))
    nbins = max(num) - min(num)
    n, bins, patches = hist(num, nbins, range=(min(num), max(num)), log=True, align='left')
    loc, scale = 0., 1.
    x = numpy.arange(bins[0], bins[-1], 1.)
    pdf = laplace(x, 0., 1.)
    plot(x, pdf)
    width = max(-min(num), max(num))
    xlim((-width, width))
    ylim((1.0, 10**7))
    show()
if __name__ == '__main__':
    main()

编辑

好的，这里尝试将其与普通拉普拉斯分布（而不是对数拉普拉斯分布）匹配。与上述尝试的不同之处：

• 直方图已经归一化
• 直方图是线性的（而不是对数的）
• laplace函数严格按照维基百科文章中指定的方式定义

输出：

如您所见，它并不是最佳匹配，但是数字（直方图和拉普拉斯PDF）至少现在在同一个范围内。我认为对数拉普拉斯会更好地匹配。我的方法（源代码在上面）没有奏效。有人能建议一种可行的方法吗？

来源：

from pylab import *
import numpy   
def laplace(x, mu, b):
    return 1.0/(2*b) * numpy.exp(-abs(x - mu)/b)
def main():
    import sys
    num = map(int, sys.stdin.read().strip().split(' '))
    nbins = max(num) - min(num)
    n, bins, patches = hist(num, nbins, range=(min(num), max(num)), log=False, align='left', normed=True)
    loc, scale = 0., 0.54
    x = numpy.arange(bins[0], bins[-1], 1.)
    pdf = laplace(x, loc, scale)
    plot(x, pdf)
    width = max(-min(num), max(num))
    xlim((-width, width))
        show()
if __name__ == '__main__':
    main()

- mpenkov

Answer 1

我已经找到了解决问题的方法。不再使用matplotlib.hist，而是在两个独立的步骤中使用numpy.histogram和matplotlib.bar来计算并绘制直方图。

我不确定是否有一种使用matplotlib.hist的方法--虽然这肯定更方便。

您可以看到它更匹配了。

我的问题现在是我需要估计PDF的scale参数。

来源:

from pylab import *
import numpy

def laplace(x, mu, b):
    """http://en.wikipedia.org/wiki/Laplace_distribution"""
    return 1.0/(2*b) * numpy.exp(-abs(x - mu)/b)

def main():
    import sys
    num = map(int, sys.stdin.read().strip().split(' '))
    nbins = max(num) - min(num)
    count, bins = numpy.histogram(num, nbins)
    bins = bins[:-1]
    assert len(bins) == nbins
    #
    # FIRST we take the log of the histogram, THEN we normalize it.
    # Clean up after divide by zero
    #
    count = numpy.log(count)
    for i in range(nbins):
        if count[i] == -numpy.inf:
            count[i] = 0
    count = count/max(count)

    loc = 0.
    scale = 4.
    x = numpy.arange(bins[0], bins[-1], 1.)
    pdf = laplace(x, loc, scale)
    pdf = pdf/max(pdf)

    width=1.0
    bar(bins-width/2, count, width=width)
    plot(x, pdf, color='r')
    xlim(min(num), max(num))
    show()

if __name__ == '__main__':
    main()

- mpenkov

---

你已经快完成了 :-). 请看我的修改建议，以便更加贴切。 - ev-br

Answer 2

1. 你的laplace()函数似乎不是拉普拉斯分布。此外，numpy.log()是自然对数（以e为底），而不是十进制。

2. 你的直方图似乎没有被归一化，而分布是被归一化的。

编辑：

1. 不要使用全局导入from pyplot import *，这会让你后悔。

2. 如果你想检查拉普拉斯分布（或其对数）的一致性，请利用后者围绕mu对称的事实：将mu固定在直方图的最大值处，你就得到了一个单参数问题。你也只能使用直方图的一半。

3. 使用numpy的直方图函数——这样你可以得到直方图本身，然后用拉普拉斯分布（和/或其对数）拟合它。卡方检验将告诉你一致性有多好（或多差）。对于拟合，你可以使用例如scipy.optimize.leastsq例程（http://www.scipy.org/Cookbook/FittingData）。

- ev-br

---

@Zhenya，感谢您的评论。为什么您说我的laplace()函数不是拉普拉斯分布？如果您查看维基百科页面，您会发现我已经完全按照定义实现了拉普拉斯概率密度函数（除了之后我将以10为底数取指数，如我原始帖子中所提到的）。关于对数底数的问题，您似乎是正确的。我被y轴上的刻度线是10的幂次方所迷惑了。它似乎确实是以e为底数。关于您的第二个观点，我会研究一下规范化直方图，谢谢。 - mpenkov

@misha：我就是不理解你为什么要用10的拉普拉斯分布幂次；也许这与你的原始数据实际有关。 - ev-br

@Zhenya：我这样做是因为我的人口概率密度函数的对数看起来像一个拉普拉斯分布（而不是通常的概率密度函数）。你看y轴的刻度是对数尺度吗？一个正常的拉普拉斯分布在这里行不通。另外，那应该真正是exp(laplace_distribution)，因为对数是以e为底的。 - mpenkov

@Zhenya：我更新了帖子，尝试匹配正常的拉普拉斯分布。也许现在我的意图更清楚了。 - mpenkov

@misha：看起来你在两张图片中使用了不同的数据集;-)。我已经更新了答案，加入了一些（几乎是随意的）想法。 - ev-br

@Zhenya：不是的，它们是同一个数据集。只是在其中一张图中y轴是线性的，在另一张图中是对数的。我想出了一个解决方法——由于原始问题变得很长，我已经将其包含在单独的答案中。你有什么意见？ - mpenkov