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

Question

Я использую 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

Answer 1

Из numpy.полифит документация, это соответствует линейной регрессии.В частности, 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

Answer 2

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

scipy.статистика.linregress

т. е.

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

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

Answer 3

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

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

Answer 4

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

Answer 5

Первоначально я опубликовал приведенные ниже тесты с целью рекомендовать 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 notebook.Объект результатов предоставляет доступ ко многим полезным статистическим показателям в дополнение к rsquared.

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

---

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

В исправительный функция, используемая в Вопросе, вычисляет коэффициент корреляции, 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 мс на цикл
• 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 раз быстрее, чем чистый python.Я не очень хорошо разбираюсь в таких вещах, как numba и pypy, поэтому кому-то другому пришлось бы заполнить эти пробелы, но я думаю, что это достаточно убедительно для меня, что corrcoef является лучшим инструментом для расчета r для простой линейной регрессии.

Вот мой код бенчмаркинга.Я скопировал вставку из записной книжки Jupyter (трудно не назвать это записной книжкой 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)
    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 6

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

Answer 7

Вот функция для вычисления взвешенный r-squared с помощью 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 - прогнозируемое значение из fit, y_{av} - среднее значение наблюдаемых данных, y_i - наблюдаемое значение данных.w_i - это взвешивание, применяемое к каждой точке данных, обычно w_i=1.SSE - это сумма квадратов из-за ошибки, а SST - общая сумма квадратов.

---

Если интересно, код в R: https://gist.github.com/dhimmel/588d64a73fa4fef02c8f (зеркало)

Answer 8

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

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

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

\сумма_{i}(y_{i} - y_bar)^2

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

Затем вы вычисляете "регрессионную сумму квадратов", которая показывает, насколько ваши ПОДОГНАННЫЕ значения отличаются от среднего

\сумма_{i}(yHat_{i} - y_bar)^2

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

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

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