Как рассчитать r-квадрат с помощью Python и Numpy?

Я использую Python и Numpy для вычисления наиболее подходящего полинома произвольной степени. Я передаю список значений x, значений y и степени полинома, который мне нужен (линейный, квадратичный и т. Д.).

Это много работает, но я также хочу вычислить r (коэффициент корреляции) и r-квадрат (коэффициент детерминации). Я сравниваю свои результаты с наиболее подходящей линией тренда в Excel и вычисляемым значением r-квадрат. Используя это, я знаю, что правильно вычисляю r-квадрат для линейного наилучшего соответствия (степень равна 1). Однако моя функция не работает для многочленов со степенью больше 1.

Excel может это сделать. Как рассчитать r-квадрат для полиномов более высокого порядка с помощью Numpy?

Вот моя функция:

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

Из документации numpy.polyfit это соответствует линейной регрессии. В частности, numpy.polyfit со степенью 'd' соответствует линейной регрессии со средней функцией

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

Очень поздний ответ, но на всякий случай кому-то для этого понадобится готовая функция:

scipy.stats.linregress

т.е.

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

как в ответе @Adam Marples.

Из yanl (еще одна-библиотека) sklearn.metricsесть r2_scoreфункция;

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

Я успешно использовал это, где 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

Первоначально я опубликовал тесты ниже с целью рекомендовать numpy.corrcoef, по глупости не понимая, что исходный вопрос уже использует corrcoefи фактически спрашивал о полиномиальных соответствиях более высокого порядка. Я добавил реальное решение вопроса о полиномиальном r-квадрате, используя statsmodels, и оставил исходные тесты, которые, хотя и не по теме, могут быть кому-то полезны.

statsmodelsимеет возможность вычислять r^2полиномиальную аппроксимацию напрямую, вот 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, следует также взглянуть на сводную информацию о подогнанной модели, которую можно распечатать или отобразить как расширенную HTML-таблицу в блокноте Jupyter / IPython. Объект результатов обеспечивает доступ ко многим полезным статистическим показателям в дополнение к rsquared.

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

Ниже мой оригинальный ответ, в котором я сравнивал различные методы линейной регрессии r ^ 2 ...

Функция corrcoef, используемая в вопросе, вычисляет коэффициент корреляции rтолько для одной линейной регрессии, поэтому она не решает вопрос о r^2подгонках полиномов более высокого порядка. Однако, как бы то ни было, я пришел к выводу, что для линейной регрессии это действительно самый быстрый и прямой метод вычисления r.

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

Это были мои результаты timeit сравнения нескольких методов для 1000 случайных (x, y) точек:

• Чистый Python (прямой rрасчет)
  • 1000 циклов, лучшее из 3: 1,59 мс на цикл
• Numpy polyfit (применимо к полиномам n-й степени)
  • 1000 петель, лучшее из 3: 326 мкс на петлю
• Руководство Numpy (прямой rрасчет)
  • 10000 циклов, лучшее из 3: 62,1 мкс на цикл
• Numpy corrcoef (прямой rрасчет)
  • 10000 петель, лучшее из 3: 56,6 мкс на петлю
• Scipy (линейная регрессия с rвыходом)
  • 1000 петель, лучшее из 3: 676 мкс на петлю
• Statsmodels (может выполнять многочлен n-й степени и многие другие варианты)
  • 1000 петель, лучшее из 3: 422 мкс на петлю

Метод corrcoef значительно превосходит вычисление r ^ 2 «вручную» с использованием методов numpy. Это более чем в 5 раз быстрее, чем метод polyfit, и примерно в 12 раз быстрее, чем scipy.linregress. Просто чтобы укрепить то, что numpy делает для вас, он в 28 раз быстрее, чем чистый питон. Я не очень разбираюсь в таких вещах, как numba и pypy, поэтому кому-то другому придется заполнить эти пробелы, но я думаю, что это достаточно убедительно для меня, что corrcoefэто лучший инструмент для расчетаr простой линейной регрессии.

Вот мой тестовый код. Я скопировал из Jupyter Notebook (трудно не назвать его IPython Notebook ...), поэтому прошу прощения, если что-то сломалось в пути. Для волшебной команды% 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)

R-квадрат - это статистика, которая применяется только к линейной регрессии.

По сути, он измеряет, насколько вариативность ваших данных может быть объяснена линейной регрессией.

Итак, вы вычисляете «общую сумму квадратов», которая представляет собой общий квадрат отклонения каждой из ваших конечных переменных от их среднего значения. . .

\ sum_ {i} (y_ {i} - y_bar) ^ 2

где y_bar - среднее значение y.

Затем вы вычисляете "сумму квадратов регрессии", которая показывает, насколько ваши FITTED значения отличаются от среднего.

\ sum_ {i} (yHat_ {i} - y_bar) ^ 2

и найдите соотношение этих двух.

Теперь все, что вам нужно сделать для полиномиального подбора, - это вставить y_hat из этой модели, но называть это r-квадрат некорректно.

Вот ссылка, которую я нашел, что немного говорит об этом.

Статья в Википедии о r-квадрате предполагает, что его можно использовать для подгонки общей модели, а не только для линейной регрессии.

Вот функция для вычисления взвешенного r-квадрата с помощью Python и Numpy (большая часть кода взята из 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 ( зеркало )

Вот очень простая функция Python для вычисления R ^ 2 из фактических и прогнозируемых значений, предполагая, что y и y_hat являются сериями панд:

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)

Из источника scipy.stats.linregress. Они используют метод средней суммы квадратов.

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

Вы можете выполнить этот код напрямую, он найдет вам полином и найдет вам 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)