如何在scipy中计算曲线拟合的可能性？

Question

我有一个非线性模型的拟合结果，长得像这样：

深色的实线是模型的拟合结果，而灰色的部分是原始数据。

简单来说，我想知道如何计算这个模型拟合的可能性，这样我就可以进行对数似然比检验。假设残差是正态分布的。

我对统计学还比较陌生，目前的想法是：

1. 从曲线拟合中得到残差，然后计算残差的方差；

2. 使用这个公式 把残差的方差代入sigma平方，x_i作为实验数据，mu作为模型拟合的结果；

3. 计算对数似然比。

有没有人能帮我解答这两个完整的问题？

1. 我的方法正确吗？（我觉得是，但确认一下会更好！）

2. 在python/scipy/statsmodels中有没有现成的函数可以帮我完成这个？

Answer 1

你的公式看起来是对的。它应该能给你和 scipy.stats.norm.logpdf(x, loc=mu, scale=sigma) 一样的结果。

既然你已经有了 mu 和 sigma 的估计值，我觉得没有一个函数可以直接把你的结果放进去进行似然比检验。

如果你有两个模型的估计值，其中一个模型是嵌套在另一个模型里的话，你可以很简单地自己计算。

http://en.wikipedia.org/wiki/Likelihood-ratio_test

这里有一个 statsmodels 中的方法部分，它计算了比较两个嵌套线性模型的 LR 检验：https://github.com/statsmodels/statsmodels/blob/master/statsmodels/regression/linear_model.py#L1531

Answer 2

你的似然函数

其实就是高斯分布的概率密度函数的对数之和。

这代表的是为你的残差拟合一个均值和标准差的可能性，而不是在给定数据的情况下你的模型的可能性。简单来说，你的方法是错误的。

因为你在做非线性最小二乘法，按照@usethedeathstar提到的，你应该直接使用F检验。考虑以下例子，修改自http://www.walkingrandomly.com/?p=5254，我们将使用R进行F检验。最后我们会讨论如何把它转到python中。

# construct the data vectors using c()
> xdata = c(-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9)
> ydata = c(0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001)
# some starting values
> p1 = 1
> p2 = 0.2
> p3 = 0.01

# do the fit
> fit1 = nls(ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata), start=list(p1=p1,p2=p2))
> fit2 = nls(ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata)+p3*xdata, start=list(p1=p1,p2=p2,p3=p3))

# summarise
> summary(fit1)

Formula: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
p1 1.881851   0.027430   68.61 2.27e-12 ***
p2 0.700230   0.009153   76.51 9.50e-13 ***
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1

Residual standard error: 0.08202 on 8 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 2.189e-06

> summary(fit2)

Formula: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata) + p3 * xdata

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
p1  1.90108    0.03520  54.002 1.96e-10 ***
p2  0.70657    0.01167  60.528 8.82e-11 ***
p3  0.02029    0.02166   0.937     0.38    
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1

Residual standard error: 0.08243 on 7 degrees of freedom

Number of iterations to convergence: 9 
Achieved convergence tolerance: 2.476e-06

> anova(fit2, fit1)
Analysis of Variance Table

Model 1: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata) + p3 * xdata
Model 2: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata)
  Res.Df Res.Sum Sq Df     Sum Sq F value Pr(>F)
1      7   0.047565                             
2      8   0.053813 -1 -0.0062473  0.9194 0.3696

这里我们有两个模型，fit1有2个参数，因此残差有8个自由度；fit2多了一个参数，残差有7个自由度。模型2真的更好吗？不，F值是0.9194，在(1,7)的自由度下，这个结果并不显著。

要得到ANOVA表：残差的自由度很简单。残差平方和：0.08202*0.08202*8=0.05381 和 0.08243*0.08243*7=0.04756293（注意：'残差标准误差：0.08243，7个自由度'，等等）。在python中，你可以通过(y_observed-y_fitted)**2来计算，因为scipy.optimize.curve_fit()不会返回残差。

F比率是0.0062473/0.047565*7，要得到P值：1-scipy.stats.f.cdf(0.9194, 1, 7)。

把它们放在一起，我们有python的等价代码：

In [1]:

import scipy.optimize as so
import scipy.stats as ss
xdata = np.array([-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9])
ydata = np.array([0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001])
def model0(x,p1,p2):
    return p1*np.cos(p2*x) + p2*np.sin(p1*x)
def model1(x,p1,p2,p3):
    return p1*np.cos(p2*x) + p2*np.sin(p1*x)+p3*x
p1, p2, p3 = 1, 0.2, 0.01
fit0=so.curve_fit(model0, xdata, ydata, p0=(p1,p2))[0]
fit1=so.curve_fit(model1, xdata, ydata, p0=(p1,p2,p3))[0]
yfit0=model0(xdata, fit0[0], fit0[1])
yfit1=model1(xdata, fit1[0], fit1[1], fit1[2])
ssq0=((yfit0-ydata)**2).sum()
ssq1=((yfit1-ydata)**2).sum()
df=len(xdata)-3
f_ratio=(ssq0-ssq1)/(ssq1/df)
p=1-ss.f.cdf(f_ratio, 1, df)
In [2]:

print f_ratio, p
0.919387419515 0.369574503394

正如@usethedeathstar指出的：当你的残差是正态分布时，非线性最小二乘法就是最大似然估计。因此F检验和似然比检验是等价的。因为F比率是似然比λ的单调变换。

或者用一种描述性的方式来看：http://www.stata.com/support/faqs/statistics/chi-squared-and-f-distributions/