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

Ответ 1

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

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

Итак, вам просто нужно вычислить R-квадрат для этой подгонки. На странице wikipedia на линейная регрессия дается полная информация. Вас интересует R ^ 2, который вы можете рассчитать несколькими способами, возможно, что easisest

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

Ответ 2

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

scipy.stats.linregress

т.е.

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

как в ответе Адама Марпла.

Ответ 3

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

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

Ответ 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

Ответ 5

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

statsmodels имеет возможность рассчитать r^2 полиномиального statsmodels, вот 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, следует также взглянуть на 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

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

• Pure Python (прямой расчет r)
  • 1000 циклов, лучшее из 3: 1,59 мс на цикл
• Полифит Numpy (применим к полиномам n-й степени)
  • 1000 петель, лучшее из 3: 326 мкс на петлю
• Numpy Manual (прямой расчет r)
  • 10000 петель, лучшее из 3: 62,1 мкс на петлю
• Numpy corrcoef (прямой расчет r)
  • 10000 петель, лучшее из 3: 56,6 мкс на петлю
• Сципи (линейная регрессия с 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 Notebook...), поэтому я прошу прощения, если что-то сломалось на пути. Команде% timeit magic требуется 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)

Ответ 6

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

Ответ 7

Вот функция для вычисления взвешенного 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

Это соответствует формула (mirror):

где f_i - предсказанное значение из подгонки, y_ {av} - среднее значение наблюдаемых данных y_i - наблюдаемое значение данных. w_i - это взвешивание, применяемое к каждой точке данных, обычно w_i = 1. SSE - сумма квадратов из-за ошибки, а SST - общая сумма квадратов.

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

Ответ 8

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

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

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

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

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

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

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

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

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

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