python - 将分布拟合到数据 : how to penalize "bad" parameter estimates?

Question

我正在使用 scipy's least-squares optimization适合 exponentially-modified gaussian distribution到一组 react 时间测量。一般来说，它运行良好，但有时，优化会偏离轨道并为参数选择一个疯狂的值——结果图显然不能很好地拟合数据。一般来说，问题似乎是由浮点精度错误引起的——我们转向 0、inf 或 nan-land。

我正在考虑做两件事:

• 使用参数同时对数据进行 CDF 和 PDF 拟合；我有两个公式。 (我正在使用核密度估计来近似数据中的 PDF。)
• 以某种方式考虑与初始参数估计值的距离(通过维基百科页面上的矩方法生成)。这些估计远非完美，但相当不错，并且似乎避开了“ float 爆炸”问题。

结合 PDF 和 CDF 拟合听起来非常简单；误差的范围甚至大体相同。但是让初始参数适合那里:我不太确定这是否是个好主意——但如果是:

• 我会如何处理规模差异？我应该将参数“错误”标准化为百分比错误吗？
• 是否有合理的方法来确定数据估计误差和参数“误差”之间的相对权重？

这些甚至是正确的问题吗？是否存在普遍认为的“正确”答案，或者“尝试一些东西直到找到似乎有效的东西”是一种好方法吗？

一个示例数据集

根据要求，这里有一个数据集，此过程对其运行不佳。我知道只有几个样本，而且数据不能很好地符合分布；我仍然希望通过优化获得“合理的外观”结果。

array([ 450.,  560.,  692.,  730.,  758.,  723.,  486.,  596.,  716.,
        695.,  757.,  522.,  535.,  419.,  478.,  666.,  637.,  569.,
        859.,  883.,  551.,  652.,  378.,  801.,  718.,  479.,  544.])

MLE 更新

在让我的 MLE 估计收敛到一个“合理”值之前，我遇到了一堆问题，直到我发现这一点:如果 X 至少包含一个 nan，np.sum(X) == nan 当 X 是一个 numpy 数组但 不是 当 X 是一个 pandas 系列时。因此，当参数开始越界时，对数似然之和正在做愚蠢的事情。

添加了一个 np.asarray() 调用，一切都很好!

Answer 1

这应该是一条评论，但我的空间用完了。

我认为最大似然拟合可能是此处最合适的方法。 ML 方法已经在 scipy.stats 中为许多发行版实现。例如，您可以调用 scipy.stats.norm.fit 找到正态分布的 MLE，并以类似的方式找到指数分布的 MLE。结合这两个生成的 MLE 参数应该可以为 Ex-Gaussian ML 拟合提供一个非常好的起始参数。事实上，我认为你的大部分数据都非常好地呈正态分布。如果是这种情况，仅正态分布的 ML 参数估计值就应该为您提供一个非常好的起始参数。

由于 Ex-Gaussian 只有 3 个参数，我认为 ML 拟合一点也不难。如果您可以提供一个数据集，而您当前的方法不能很好地处理该数据集，那么展示一个真实的例子会更容易。

好的，给你:

>>> import scipy.special as sse
>>> import scipy.stats as sss
>>> import scipy.optimize as so
>>> from numpy import *

>>> def eg_pdf(p, x): #defines the PDF
    m=p[0]
    s=p[1]
    l=p[2]
    return 0.5*l*exp(0.5*l*(2*m+l*s*s-2*x))*sse.erfc((m+l*s*s-x)/(sqrt(2)*s))

>>> xo=array([ 450.,  560.,  692.,  730.,  758.,  723.,  486.,  596.,  716.,
        695.,  757.,  522.,  535.,  419.,  478.,  666.,  637.,  569.,
        859.,  883.,  551.,  652.,  378.,  801.,  718.,  479.,  544.])

>>> sss.norm.fit(xo) #get the starting parameter vector form the normal MLE
(624.22222222222217, 132.23977474531389)

>>> def llh(p, f, x): #defines the negative log-likelihood function
    return -sum(log(f(p,x)))

>>> so.fmin(llh, array([624.22222222222217, 132.23977474531389, 1e-6]), (eg_pdf, xo)) #yeah, the data is not good
Warning: Maximum number of function evaluations has been exceeded.
array([  6.14003407e+02,   1.31843250e+02,   9.79425845e-02])

>>> przt=so.fmin(llh, array([624.22222222222217, 132.23977474531389, 1e-6]), (eg_pdf, xo), maxfun=1000) #so, we increase the number of function call uplimit
Optimization terminated successfully.
         Current function value: 170.195924
         Iterations: 376
         Function evaluations: 681

>>> llh(array([624.22222222222217, 132.23977474531389, 1e-6]), eg_pdf, xo)
400.02921290185645
>>> llh(przt, eg_pdf, xo) #quite an improvement over the initial guess
170.19592431051217
>>> przt
array([  6.14007039e+02,   1.31844654e+02,   9.78934519e-02])

此处使用的优化器(fmin，或 Nelder-Mead 单纯形算法)不使用任何来自梯度的信息，并且通常比使用的优化器运行得慢得多。看来指数高斯的负对数似然函数的导数可以很容易地写成一个封闭的形式。如果是这样，利用梯度/导数的优化器将是更好、更有效的选择(例如 fmin_bfgs)。

要考虑的另一件事是参数约束。根据定义，对于指数高斯分布，sigma 和 lambda 必须为正。您可以使用约束优化器(例如 fmin_l_bfgs_b)。或者，您可以优化:

>>> def eg_pdf2(p, x): #defines the PDF
    m=p[0]
    s=exp(p[1])
    l=exp(p[2])
    return 0.5*l*exp(0.5*l*(2*m+l*s*s-2*x))*sse.erfc((m+l*s*s-x)/(sqrt(2)*s))

由于MLE的函数不变性，这个函数的MLE应该和原来的eg_pdf一样。除了 exp() 之外，您还可以使用其他转换将 (-inf, +inf) 投影到 (0, +inf).

你也可以考虑http://en.wikipedia.org/wiki/Lagrange_multiplier .