如何使用scipy.optimize.curve_fit()函数返回的结果计算决定系数R2

问题引入

当我们需要对一批数据做曲线拟合的时候，来自python的scipy包下的curve_fit()函数往往是一个不错的选择，但curve_fit()函数返回的结果只有拟合曲线的参数popt和参数的估计协方差pcov（etismatated covarianve of popt）[1]。而作为回归模型重要的衡量参数——决定系数\(R^2\)，并没有作为函数的返回结果直接提供。（函数的具体使用参见官方doc）

接下来将解说如何利用curve_fit()函数返回的结果来计算\(R^2\)。

背景知识

决定系数\(R^2\)的定义

先放一段维基百科中的定义[2]

In statistics, the coefficient of determination, denoted \(R^2\) or \(r^2\) and pronounced "R squared", is the proportion of the variation in the dependent variable that is predictable from the independent variable(s).

翻译到中文就是：

在统计学中，决定系数（记作\(R^2\)或\(r^2\)）是用于度量因变量的差异中可由自变量解释的部分所占的比例。

我们一般使用决定系数来判断回归模型的解释力，通俗一点说就是决定系数越高的模型，拟合度越高。

决定系数\(R^2\)的计算

定义一个数据集包括\(y_1,\cdots,y_n\)共\(n\)个观测值，相应的经过模型拟合后的结果\(f_1,\cdots,f_n\)，同时定义残差为\(e_i=y_i-f_i\)。[2]

若\(\overline{y}\)为观测值的均值，定义如下：

\[\overline{y}=\frac{1}{n}\sum_{i=1}^{n}{y_i} \tag{1} \]

于是（拟合结果与）数据集的差异可以用以下两个平方和公式来衡量：

• 总体平方和（与数据的方差成正比）

  \[SS_{tot}=\sum_{i}{(y_i-\overline{y})}^{2} \tag{2} \]

• 残差平方和

  \[SS_{res}=\sum_{i}{(y_i-f_i)}^{2}=\sum_{i}{e_i^2} \tag{3} \]

由此，决定系数的最常用的定义如下：

\[R^2=1-\frac{SS_{res}}{SS_{tot}} \tag{4} \]

在理想状态下，模型拟合的结果和观测值正好完全一致，这样就能计算得到\(SS_{res}=0\)和\(R^2=0\)。作为基准，如果模型一直拟合出结果为\(\overline{y}\)，那么其\(R^2=0\)。具有比此基准更差的预测结果的模型，其\(R^2\)将为负值。

解决问题

示例数据

首先定义一个待拟合的数据集。

在\(y=x\)的基础上，给每个\(y\)随机加减\([-2,2]\)内的任意值，由此构建出示例数据，代码如下：

import numpy as np

# define func: y=ax+b
def func(x, a, b):
    x = np.array(x)
    return a * x + b

X = range(0, 50) # x
y = func(X, a=1, b=0) # baseline function: y=x
res = np.random.uniform(low=-2.0, high=2.0, size=len(X))
y = [y[i] + res[i] for i in range(0, len(y))]

绘制出散点图如下：

使用curve_fit()函数做拟合

在示例数据和之前定义的函数基础上，我们对数据做拟合，代码实现如下：

from scipy.optimize import curve_fit

popt, pcov = curve_fit(func, X, y)

拟合出的函数图像如下：

根据公式计算拟合曲线的决定系数

依照前面的公式，我们可以将决定系数的计算分为以下几步：

1. 计算均值
2. 计算总体平方和
3. 计算残差平方和
4. 计算决定系数

代码实现如下[3]：

import numpy as np

mean = np.mean(y)  # 1.y mean
ss_tot = np.sum((y - mean) ** 2)  # 2.total sum of squares
ss_res = np.sum((y - func(X, *popt)) ** 2)  # 3.residual sum of squares
r_squared = 1 - (ss_res / ss_tot)  # 4.r squared

我们将前面的结果和计算出的决定系数绘制到同一张图中：

参考

[1] scipy.optimize.curve_fit — SciPy v1.6.3 Reference Guide

[2] Coefficient of determination - Wikipedia

[3] python - Getting the r-squared value using curve_fit - Stack Overflow