拟合某些参数固定的双峰高斯分布

Question

问题： 我想将经验数据拟合到双峰正态分布，从物理环境中我知道峰的距离（固定）以及两个峰必须具有相同的标准差。

我试图用 scipy.stats.rv_continous 创建一个自己的分布（见下面的代码），但参数总是适合 1。有人知道发生了什么，或者可以给我指出一个不同的方法来解决问题了吗？

详情： 我避免了 loc 和 scale 参数并将它们实现为 m 和 s 直接进入_pdf-方法因为峰距delta不受scale的影响。为了弥补这一点，我在 fit 方法中将它们固定为 floc=0 和 fscale=1 并且实际上想要 m、s 和权重的拟合参数山峰 w

我在示例数据中期望的是峰值分布在 x=-450 和 x=450 (=> m=0) 附近。 stdev s 应该在 100 或 200 左右，但不是 1.0，权重 w 应该是大约。 0.5

from __future__ import division
from scipy.stats import rv_continuous
import numpy as np


class norm2_gen(rv_continuous):
    def _argcheck(self, *args):
        return True

    def _pdf(self, x, m, s, w, delta):
        return np.exp(-(x-m+delta/2)**2 / (2. * s**2)) / np.sqrt(2. * np.pi * s**2) * w + \
               np.exp(-(x-m-delta/2)**2 / (2. * s**2)) / np.sqrt(2. * np.pi * s**2) * (1 - w)


norm2 = norm2_gen(name='norm2')

data = [487.0, -325.5, -159.0, 326.5, 538.0, 552.0, 563.0, -156.0, 545.5, 341.0, 530.0, -156.0, 473.0, 328.0, -319.5, -287.0, -294.5, 153.5, -512.0, 386.0, -129.0, -432.5, -382.0, -346.5, 349.0, 391.0, 299.0, 364.0, -283.0, 562.5, -42.0, 214.0, -389.0, 42.5, 259.5, -302.5, 330.5, -338.0, 508.5, 319.5, -356.5, 421.5, 543.0]

m, s, w, delta, loc, scale = norm2.fit(data, fdelta=900, floc=0, fscale=1)
print m, s, w, delta, loc, scale

>>> 1.0 1.0 1.0 900 0 1

Answer 1

经过一些调整后，我能够使您的分布适合数据：

像你一样使用w，你有一个隐含的约束0 <= w <= 1。fit()方法使用的求解器不知道这一点约束，因此 w 可能会得到不合理的值。处理此类约束的一种方法是允许 w 为任意实数值，但在 PDF 的公式中，将 w 转换为 0 和 1 之间的分数 phi 使用phi = 0.5 + arctan(w)/pi.
通用 fit() 方法使用数值优化例程来查找最大似然估计。与大多数此类例程一样，它需要优化的起点。默认起点全为 1，但这并不总是有效。您可以通过在数据后向 fit() 提供值作为位置参数来选择不同的起点。我在脚本中使用的值有效；我没有探讨结果对这些起始值有多敏感。

我做了两个估计。在第一个中，我让 delta 成为一个自由参数，在第二个中，我将 delta 固定为 900。

下面的脚本生成以下图：

这是脚本：

from __future__ import division
from scipy.stats import rv_continuous
import numpy as np
import matplotlib.pyplot as plt


class norm2_gen(rv_continuous):
    def _argcheck(self, *args):
        return True

    def _pdf(self, x, m, s, w, delta):
        phi = 0.5 + np.arctan(w)/np.pi
        return np.exp(-(x-m+delta/2)**2 / (2. * s**2)) / np.sqrt(2. * np.pi * s**2) * phi + \
               np.exp(-(x-m-delta/2)**2 / (2. * s**2)) / np.sqrt(2. * np.pi * s**2) * (1 - phi)

norm2 = norm2_gen(name='norm2')


data = [487.0, -325.5, -159.0, 326.5, 538.0, 552.0, 563.0, -156.0, 545.5,
        341.0, 530.0, -156.0, 473.0, 328.0, -319.5, -287.0, -294.5, 153.5,
        -512.0, 386.0, -129.0, -432.5, -382.0, -346.5, 349.0, 391.0, 299.0,
        364.0, -283.0, 562.5, -42.0, 214.0, -389.0, 42.5, 259.5, -302.5,
        330.5, -338.0, 508.5, 319.5, -356.5, 421.5, 543.0]

# In the fit method, the positional arguments after data are the initial
# guesses that are passed to the optimization routine that computes the MLE.
# First let's see what we get if delta is not fixed.
m, s, w, delta, loc, scale = norm2.fit(data, 1.0, 1.0, 0.0, 900.0, floc=0, fscale=1)

# Fit the disribution with delta fixed.
fdelta = 900
m1, s1, w1, delta1, loc, scale = norm2.fit(data, 1.0, 1.0, 0.0, fdelta=fdelta, floc=0, fscale=1)

plt.hist(data, bins=12, normed=True, color='c', alpha=0.65)
q = np.linspace(-800, 800, 1000)
p = norm2.pdf(q, m, s, w, delta)
p1 = norm2.pdf(q, m1, s1, w1, fdelta)
plt.plot(q, p, 'k', linewidth=2.5, label='delta=%6.2f (fit)' % delta)
plt.plot(q, p1, 'k--', linewidth=2.5, label='delta=%6.2f (fixed)' % fdelta)
plt.legend(loc='best')
plt.show()

拟合某些参数固定的双峰高斯分布

Fitting bimodal gaussian distribution with some parameters fixed

python

statistics

distribution

scipy