使用对数拉普拉斯分布函数绘制Matplotlib直方图

Question

（在深入阅读源代码之前，请务必查看帖子末尾的编辑内容）

我正在绘制一个看起来符合对数拉普拉斯分布的人口直方图：

我想为这个直方图画一条最佳拟合线，以验证我的假设，但我在得到有意义的结果时遇到了问题。

我使用了维基百科上的拉普拉斯概率密度函数（PDF）定义，并将其结果取10的幂（这样做是为了“逆转”对数直方图的影响）。

我哪里做错了呢？

这是我的代码。我通过标准输入来处理数据（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()

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()

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）。