如何使用 Python 和 Numpy 计算 r 平方？

小开

R 平方是一个只适用于线性回归的统计量。

本质上，它衡量的是你的数据中有多少变化可以用线性回归来解释。

所以，你计算“总平方和”，这是总平方偏差的每个结果变量从他们的平均值. 。

$\sum_{i}(y_{i} - y_{bar})^2$

其中 y _ bar 是 y 的平均值。

然后，计算“回归平方和”，即 FITTED 值与平均值的差异

$\sum_{i}(yhat_{i} - y_{bar})^2$

找出这两者的比例。

现在，你需要做的就是把这个模型的 y _ hat 插入到多项式拟合中，但是把它叫做 r 平方是不准确的。

在这里是我发现的一个链接，说了一点。

小开

最佳答案

从 Numpy Polyfit的文档来看，这是一个合适的线性回归。具体来说，度数为 d 的 numpy.polyfit 符合平均值函数的线性回归

E (y | x) = p _ d * x * * d + p _ { d-1} * x * * (d-1) + ... + p _ 1 * x + p _ 0

所以你只需要计算这个拟合的 R 平方。线性回归的维基百科页面提供了详细信息。您对 R ^ 2感兴趣，可以用几种方法计算，最简单的方法可能是

SST = Sum(i=1..n) (y_i - y_bar)^2
SSReg = Sum(i=1..n) (y_ihat - y_bar)^2
Rsquared = SSReg/SST

其中我使用‘ y _ bar’表示 y 的平均值，‘ y _ ihat’表示每个点的拟合值。

我不是很熟悉 numpy (我通常使用 R) ，所以可能有一个更简洁的方法来计算 R 的平方，但是下面应该是正确的

import numpy


# Polynomial Regression
def polyfit(x, y, degree):
results = {}


coeffs = numpy.polyfit(x, y, degree)


# Polynomial Coefficients
results['polynomial'] = coeffs.tolist()


# r-squared
p = numpy.poly1d(coeffs)
# fit values, and mean
yhat = p(x)                         # or [p(z) for z in x]
ybar = numpy.sum(y)/len(y)          # or sum(y)/len(y)
ssreg = numpy.sum((yhat-ybar)**2)   # or sum([ (yihat - ybar)**2 for yihat in yhat])
sstot = numpy.sum((y - ybar)**2)    # or sum([ (yi - ybar)**2 for yi in y])
results['determination'] = ssreg / sstot


return results

小开

维基百科上关于 R 平方的文章表明，它可能用于一般模型的拟合，而不仅仅是线性回归。

小开

一个非常晚的回复，但只是为了以防有人需要一个现成的功能:

Scipy.stats.linregress

也就是说。

slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(x, y)

就像“亚当 · 马普尔斯的回答”一样。

小开

我已经成功地使用了它，其中 x 和 y 是类似于数组的。

注: 只供线性回归使用

def rsquared(x, y):
""" Return R^2 where x and y are array-like."""


slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(x, y)
return r_value**2

小开

来自 yanl (另一个库)的 sklearn.metrics具有 r2_score功能;

from sklearn.metrics import r2_score


coefficient_of_dermination = r2_score(y, p(x))

小开

我最初发布下面的基准是为了推荐 numpy.corrcoef，愚蠢地没有意识到最初的问题已经使用了 corrcoef，实际上是在询问更高阶的多项式拟合。我已经使用 statsmodel 为多项式 r 平方问题添加了一个实际的解决方案，并且留下了最初的基准，虽然离题了，但对某些人来说可能是有用的。

statsmodels具有直接计算多项式拟合的 r^2的能力，这里有两种方法。

import statsmodels.api as sm
import statsmodels.formula.api as smf


# Construct the columns for the different powers of x
def get_r2_statsmodels(x, y, k=1):
xpoly = np.column_stack([x**i for i in range(k+1)])
return sm.OLS(y, xpoly).fit().rsquared


# Use the formula API and construct a formula describing the polynomial
def get_r2_statsmodels_formula(x, y, k=1):
formula = 'y ~ 1 + ' + ' + '.join('I(x**{})'.format(i) for i in range(1, k+1))
data = {'x': x, 'y': y}
return smf.ols(formula, data).fit().rsquared # or rsquared_adj

为了进一步利用 statsmodels的优势，我们还应该看看合适的模型摘要，它可以打印或显示为 Jupyter/IPython 笔记本中的一个丰富的 HTML 表格。Result 对象除了 rsquared之外，还提供了对许多有用的统计指标的访问。

model = sm.OLS(y, xpoly)
results = model.fit()
results.summary()

下面是我的原始答案，我在其中对多种线性回归 r ^ 2方法进行了基准测试... ..。

问题中使用的 Corcoef函数只计算一个线性回归的相关系数 r，所以它并没有解决高阶多项式拟合的 r^2问题。然而，无论如何，我发现对于线性回归来说，这确实是计算 r最快、最直接的方法。

def get_r2_numpy_corrcoef(x, y):
return np.corrcoef(x, y)[0, 1]**2

这是我对1000个随机(x，y)点的一系列方法进行比较得出的时间结果:

纯 Python (直接 r计算)
- 1000个循环，最好的3:1.59毫秒每个循环
粗糙多项式拟合(适用于 n 次多项式拟合)
- 1000个循环，最好是每个循环3:326μs
麻木手册(直接 r计算)
- 10000个循环，最好是每个循环3:62.1 μs
麻木的脑袋(直接 r计算)
- 10000个循环，最好是每个循环3:56.6 μs
Scypy (线性回归输出为 r)
- 1000个循环，最好是每个循环3:676μs
统计模型(可以做 n 次多项式和许多其他拟合)
- 1000个循环，最好是每个循环3:422μs

Corrcoef 方法比使用 numpy 方法“手动”计算 r ^ 2要好得多。它比 polyfit 方法快5倍以上，比 scypy.linregress 快约12倍。只是为了强调一下 numpy 为你做了什么，它比纯 Python 快了28倍。我并不精通 numba 和 py 之类的东西，所以需要其他人来填补这些空白，但我认为这足以让我相信，对于简单线性回归来说，corrcoef是计算 r的最佳工具。

这是我的基准测试代码。我复制粘贴从木星笔记本(很难不称之为 IPython 笔记本...) ，所以我道歉，如果有什么打破的方式。% timeit 魔法命令需要 IPython。

import numpy as np
from scipy import stats
import statsmodels.api as sm
import math


n=1000
x = np.random.rand(1000)*10
x.sort()
y = 10 * x + (5+np.random.randn(1000)*10-5)


x_list = list(x)
y_list = list(y)


def get_r2_numpy(x, y):
slope, intercept = np.polyfit(x, y, 1)
r_squared = 1 - (sum((y - (slope * x + intercept))**2) / ((len(y) - 1) * np.var(y, ddof=1)))
return r_squared
    

def get_r2_scipy(x, y):
_, _, r_value, _, _ = stats.linregress(x, y)
return r_value**2
    

def get_r2_statsmodels(x, y):
return sm.OLS(y, sm.add_constant(x)).fit().rsquared
    

def get_r2_python(x_list, y_list):
n = len(x_list)
x_bar = sum(x_list)/n
y_bar = sum(y_list)/n
x_std = math.sqrt(sum([(xi-x_bar)**2 for xi in x_list])/(n-1))
y_std = math.sqrt(sum([(yi-y_bar)**2 for yi in y_list])/(n-1))
zx = [(xi-x_bar)/x_std for xi in x_list]
zy = [(yi-y_bar)/y_std for yi in y_list]
r = sum(zxi*zyi for zxi, zyi in zip(zx, zy))/(n-1)
return r**2
    

def get_r2_numpy_manual(x, y):
zx = (x-np.mean(x))/np.std(x, ddof=1)
zy = (y-np.mean(y))/np.std(y, ddof=1)
r = np.sum(zx*zy)/(len(x)-1)
return r**2
    

def get_r2_numpy_corrcoef(x, y):
return np.corrcoef(x, y)[0, 1]**2
    

print('Python')
%timeit get_r2_python(x_list, y_list)
print('Numpy polyfit')
%timeit get_r2_numpy(x, y)
print('Numpy Manual')
%timeit get_r2_numpy_manual(x, y)
print('Numpy corrcoef')
%timeit get_r2_numpy_corrcoef(x, y)
print('Scipy')
%timeit get_r2_scipy(x, y)
print('Statsmodels')
%timeit get_r2_statsmodels(x, y)

7/28/21 Benchmark result (Python 3.7，numpy 1.19，scypy 1.6，statsmodel 0.12)

Python
2.41 ms ± 180 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Numpy polyfit
318 µs ± 44.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Numpy Manual
79.3 µs ± 4.05 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
Numpy corrcoef
83.8 µs ± 1.37 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
Scipy
221 µs ± 7.12 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Statsmodels
375 µs ± 3.63 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

小开

下面是一个用 Python 和 Numpy 计算 加重了 r 平方的函数(大部分代码来自 sklearn) :

from __future__ import division
import numpy as np


def compute_r2_weighted(y_true, y_pred, weight):
sse = (weight * (y_true - y_pred) ** 2).sum(axis=0, dtype=np.float64)
tse = (weight * (y_true - np.average(
y_true, axis=0, weights=weight)) ** 2).sum(axis=0, dtype=np.float64)
r2_score = 1 - (sse / tse)
return r2_score, sse, tse

例如:

from __future__ import print_function, division
import sklearn.metrics


def compute_r2_weighted(y_true, y_pred, weight):
sse = (weight * (y_true - y_pred) ** 2).sum(axis=0, dtype=np.float64)
tse = (weight * (y_true - np.average(
y_true, axis=0, weights=weight)) ** 2).sum(axis=0, dtype=np.float64)
r2_score = 1 - (sse / tse)
return r2_score, sse, tse


def compute_r2(y_true, y_predicted):
sse = sum((y_true - y_predicted)**2)
tse = (len(y_true) - 1) * np.var(y_true, ddof=1)
r2_score = 1 - (sse / tse)
return r2_score, sse, tse


def main():
'''
Demonstrate the use of compute_r2_weighted() and checks the results against sklearn
'''
y_true = [3, -0.5, 2, 7]
y_pred = [2.5, 0.0, 2, 8]
weight = [1, 5, 1, 2]
r2_score = sklearn.metrics.r2_score(y_true, y_pred)
print('r2_score: {0}'.format(r2_score))
r2_score,_,_ = compute_r2(np.array(y_true), np.array(y_pred))
print('r2_score: {0}'.format(r2_score))
r2_score = sklearn.metrics.r2_score(y_true, y_pred,weight)
print('r2_score weighted: {0}'.format(r2_score))
r2_score,_,_ = compute_r2_weighted(np.array(y_true), np.array(y_pred), np.array(weight))
print('r2_score weighted: {0}'.format(r2_score))


if __name__ == "__main__":
main()
#cProfile.run('main()') # if you want to do some profiling

产出:

r2_score: 0.9486081370449679
r2_score: 0.9486081370449679
r2_score weighted: 0.9573170731707317
r2_score weighted: 0.9573170731707317

这对应于配方奶粉(镜子) :

在 f _ i 为拟合预测值的情况下，y _ { av }为观测数据的均值，y _ i 为观测数据的均值。W _ i 是应用于每个数据点的权重，通常是 w _ i = 1。SSE 是由于误差导致的平方和，SST 是平方的总和。

如果感兴趣，请查看 R: https://gist.github.com/dhimmel/588d64a73fa4fef02c8f(镜子)中的代码

小开

他们使用的是平均平方和法。

import numpy as np


x = np.array(x)
y = np.array(y)


# average sum of squares:
ssxm, ssxym, ssyxm, ssym = np.cov(x, y, bias=1).flat


r_num = ssxym
r_den = np.sqrt(ssxm * ssym)
r = r_num / r_den


if r_den == 0.0:
r = 0.0
else:
r = r_num / r_den


if r > 1.0:
r = 1.0
elif r < -1.0:
r = -1.0

小开

下面是一个非常简单的 python 函数，它根据假设 y 和 y _ hat 是熊猫级数的实际值和预测值计算 R ^ 2:

def r_squared(y, y_hat):
y_bar = y.mean()
ss_tot = ((y-y_bar)**2).sum()
ss_res = ((y-y_hat)**2).sum()
return 1 - (ss_res/ss_tot)

小开

你可以直接执行这段代码，它会找到多项式，并且会找到 R 值 ，如果你需要更多的解释，你可以在下面写一个注释。

from scipy.stats import linregress
import numpy as np


x = np.array([1,2,3,4,5,6])
y = np.array([2,3,5,6,7,8])


p3 = np.polyfit(x,y,3) # 3rd degree polynomial, you can change it to any degree you want
xp = np.linspace(1,6,6)  # 6 means the length of the line
poly_arr = np.polyval(p3,xp)


poly_list = [round(num, 3) for num in list(poly_arr)]
slope, intercept, r_value, p_value, std_err = linregress(x, poly_list)
print(r_value**2)

小开

使用 numpy 模块(在 python3中测试) :

import numpy as np
def linear_regression(x, y):
coefs = np.polynomial.polynomial.polyfit(x, y, 1)
ffit = np.poly1d(coefs)
m = ffit[0]
b = ffit[1]
eq = 'y = {}x + {}'.format(round(m, 3), round(b, 3))
rsquared = np.corrcoef(x, y)[0, 1]**2
return rsquared, eq, m, b


rsquared, eq, m, b = linear_regression(x,y)
print(rsquared, m, b)
print(eq)

产出:

0.013378252355751777 0.1316331351105754 0.7928782850418713
y = 0.132x + 0.793

注: R2≠ R2
R2被称为“决定系数”
R2是皮尔逊系数的平方

正式合并为 r2的 R2可能是你想要的，因为它是最小二乘拟合，比 r2的简单分数要好。Numpy 并不害怕称之为“ corrcoef”，它假设 Pearson 是事实上的相关系数。