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

Question

我正在使用Python和Numpy计算任意度数的最佳拟合多项式。我传递了x值，y值以及要拟合的多项式的阶数（线性，二次等）的列表。

这很有效，但我还想计算r（相关系数）和r-平方（确定系数）。我正在将我的结果与Excel的最佳拟合趋势线功能及其计算的r平方值进行比较。使用这个，我知道我正在为线性最佳拟合（度等于1）正确计算r平方。但是，我的函数不适用于次数大于1的多项式。

Excel可以做到这一点。如何使用Numpy计算高阶多项式的r平方？

这是我的职能：

import numpy

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

    coeffs = numpy.polyfit(x, y, degree)
     # Polynomial Coefficients
    results['polynomial'] = coeffs.tolist()

    correlation = numpy.corrcoef(x, y)[0,1]

     # r
    results['correlation'] = correlation
     # r-squared
    results['determination'] = correlation**2

    return results

Answer 1

根据numpy.polyfit文档，它适合线性回归。具体来说，度为'd'的numpy.polyfit与均值函数拟合线性回归]

E（y | x）= p_d * x ** d + p_ {d-1} * x **（d-1）+ ... + p_1 * x + p_0

因此您只需要计算该拟合的R平方即可。 linear regression上的Wikipedia页面提供了完整的详细信息。您对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

Answer 2

回复很晚，但以防万一有人需要为此准备就绪的功能：

scipy.stats.linregress

即

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

如@Adam Marples的回答。

Answer 3

从yanl（另一个图书馆）sklearn.metrics具有r2_square功能；

r2_square

Answer 4

我一直在成功使用它，其中x和y类似于数组。

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

Answer 5

我最初在下面发布基准测试是为了推荐

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

，但愚蠢地没有意识到原来的问题已经使用numpy.corrcoef，实际上是在询问高阶多项式拟合。我已经使用statsmodels为多项式r平方问题添加了实际的解决方案，并且留下了原始基准测试，虽然离题，但对某人可能很有用。

---

[corrcoef可以直接计算多项式拟合的statsmodels，这里有2种方法...

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

，还应该查看拟合的模型摘要，可以在Jupyter / IPython笔记本中将其打印或显示为丰富的HTML表。除了statsmodels，结果对象还提供对许多有用的统计指标的访问。

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

Question中使用的

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

函数仅针对单个线性回归计算相关系数corrcoef，因此对于高阶多项式拟合，它没有解决r问题。但是，对于它的价值而言，我发现对于线性回归，它确实是计算r^2的最快，最直接的方法。

r
这些是我通过比较一堆方法对1000个随机（x，y）点的时间结果：

• 纯Python（直接

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

  计算）
  • 1000个循环，每个循环最好3：1.59毫秒
  • Numpy多项式拟合（适用于n次多项式拟合）
    • 1000个循环，每个循环最好3：326 µs
• Numpy手册（直接计算r）
  • 10000个循环，最好为3：每个循环62.1 µs
• Numpy corrcoef（直接计算r）
  • 10000个循环，每个循环中最好为3：56.6 µs
• Scipy（以r作为输出的线性回归）
  • 1000个循环，每个循环最好3：676 µs
• Statsmodels（可以执行n次多项式和许多其他拟合）
  • 1000个循环，每个循环最好3：422 µs

Corrcoef方法几乎胜过使用numpy方法“手动”计算r ^ 2的方法。它比polyfit方法快5倍以上，比scipy.linregress快12倍左右。只是为了增强numpy为您所做的工作，它比纯python快28倍。我不精通numba和pypy之类的东西，因此必须由其他人来填补这些空白，但是我认为这足以说服r是计算简单线性的corrcoef的最佳工具回归。

这是我的基准测试代码。我从Jupyter笔记本复制粘贴（很难不称它为IPython笔记本...），因此如果途中发生任何问题，我深表歉意。 ％timeit magic命令需要IPython。

r

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)
    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)

上的维基百科文章建议将其用于一般模型拟合，而不仅仅是线性回归。

这里是使用Python和Numpy计算weighted

r平方的函数（大多数代码来自sklearn）：

r-squareds
示例：
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
（formula）：
mirror

其中f_i是拟合的预测值，y_ {av}是观测数据的平均值y_i是观测数据值。 w_i是应用于每个数据点的权重，通常w_i = 1。 SSE是由于误差引起的平方和，而SST是平方和。

---

如果有兴趣，R中的代码：（https://gist.github.com/dhimmel/588d64a73fa4fef02c8f）

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

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

因此，您将计算“平方和”，它是每个结果变量与其平均值的总平方偏差。 。 。

\ sum_ {i}（y_ {i}-y_bar）^ 2

其中y_bar是y的平均值。

然后，您计算出“平方的平方和”，即您的FITTED值与均值相差多少

\ sum_ {i}（yHat_ {i}-y_bar）^ 2

并找到这两个比率。

现在，要进行多项式拟合，您所要做的就是插入该模型的y_hat，但称该r平方并不准确。

mirror是我发现的一个链接，它对此有所说明。

来自scipy.stats.linregress源。他们使用平均平方和方法。

Here

Answer 6

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)
    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)

上的维基百科文章建议将其用于一般模型拟合，而不仅仅是线性回归。

Answer 7

这里是使用Python和Numpy计算weighted

Answer 8

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

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

Answer 9

来自scipy.stats.linregress源。他们使用平均平方和方法。

Here