将概率分布拟合到数据并找到其累积分布函数

Question

我想把一种不对称的概率分布适配到我的数据上，我觉得指数修正的高斯分布可能是个不错的选择。

m=array([ 16.25,  16.75,  17.25,  17.75,  18.25,  18.75,  19.25,  19.75,
          20.25,  20.75,  21.25,  21.75,  22.25,  22.75,  23.25,  23.75,
          24.25,  24.75,  25.25,  25.75,  26.25,  26.75,  27.25,  27.75,
          28.25,  28.75,  29.25,  29.75,  30.25,  30.75])
pdf=array([  0.00000000e+00,   2.40818784e-04,   1.38470801e-03,
             1.62552679e-03,   3.07043949e-03,   3.37146297e-03,
             5.47862733e-03,   8.36845274e-03,   1.61348585e-02,
             1.92052980e-02,   2.79951836e-02,   3.97953040e-02,
             4.95484648e-02,   7.09211318e-02,   9.50030102e-02,
             1.40878989e-01,   1.90186635e-01,   2.42022878e-01,
             2.77302830e-01,   2.69054786e-01,   2.40397351e-01,
             1.74593618e-01,   9.16917520e-02,   2.41420831e-02,
             7.22456352e-03,   3.01023480e-04,   0.00000000e+00,
             0.00000000e+00,   0.00000000e+00,   6.02046960e-05])

我想使用scipy.optimize这个库，同时可以控制拟合的好坏，看看能不能通过改变输入参数的初始条件来改善卡方值。我写了以下这段代码：

import scipy.special as sse
from math import *
import numpy as np
import scipy.optimize
#defines the PDF of an exponentially modified Gaussian distribution
fitfunc =lambda p,x: 0.5*p[2]*np.exp(0.5*p[2]*(2*p[0]+p[2]*p[1]*p[1]-2*x))*sse.erfc((p[0]+p[2]*p[1]*p[1]-x)/(np.sqrt(2)*p[1]))
"""Deviations of data from fitted curve"""
errfunc = lambda p, x, y: fitfunc(p, x) - y
#initial values
p0=[24,1,1]
p1, success = scipy.optimize.leastsq(errfunc, p0, args=(pdf, m), maxfev=10000)

更新: 嗯，我刚刚选择了numpy.exp，第一个问题解决了，但leastsq还是给我不可靠的输出，我该怎么办呢？另外，我还想得到这个分布的CDF。

Answer 1

最小二乘法并不是用来将数据拟合到特定概率分布函数（pdf）上的正确方法。

你可能想要使用的是最大似然法，也就是说，你要最大化 p(x | a)，其中 a 是分布的参数，x 是数据。通常，我们会先计算对数似然，并假设数据是独立的，所以 log p(x | a) = sum(log(pdf(x[i], a)), i, 1, n)。

你需要使用一个最小化函数，把 log p(x | a) 作为需要最小化的函数，而 a 则是它的自由参数。