Jak obliczyć r-kwadrat za pomocą Pythona i Numpy?

Używam Pythona i Numpy, aby obliczyć najlepiej dopasowany wielomian dowolnego stopnia. Przekazuję listę wartości x, wartości y i stopień wielomianu, który chcę dopasować (liniowy, kwadratowy itp.).

To dużo działa, ale chcę też obliczyć r (współczynnik korelacji) i r-kwadrat (współczynnik determinacji). Porównuję swoje wyniki z funkcją najlepiej dopasowanej linii trendu programu Excel i obliczaną przez nią wartością r-kwadrat. Używając tego, wiem, że poprawnie obliczam r-kwadrat dla najlepszego dopasowania liniowego (stopień równy 1). Jednak moja funkcja nie działa dla wielomianów o stopniu większym niż 1.

Excel to potrafi. Jak obliczyć r-kwadrat dla wielomianów wyższego rzędu za pomocą Numpy?

Oto moja funkcja:

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

Z dokumentacji numpy.polyfit pasuje regresja liniowa. Konkretnie, numpy.polyfit z stopniem „d” pasuje do regresji liniowej z funkcją średniej

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

Wystarczy więc obliczyć R-kwadrat dla tego dopasowania. Szczegółowe informacje można znaleźć na stronie wikipedii poświęconej regresji liniowej . Interesuje Cię R ^ 2, który możesz obliczyć na kilka sposobów, prawdopodobnie najłatwiejszy

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

Gdzie używam „y_bar” jako średniej y, a „y_ihat” jako wartości dopasowania dla każdego punktu.

Nie jestem zbyt zaznajomiony z numpy (zwykle pracuję w R), więc prawdopodobnie istnieje bardziej uporządkowany sposób obliczenia R-kwadrat, ale poniższe powinny być poprawne

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

Bardzo późna odpowiedź, ale na wypadek, gdyby ktoś potrzebował do tego gotowej funkcji:

scipy.stats.linregress

to znaczy

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

jak w odpowiedzi @Adam Marples.

From yanl (yet-another-library) sklearn.metricsma r2_scorefunkcję;

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

Z powodzeniem używam tego, gdzie x i y są podobne do tablicy.

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

Pierwotnie zamieściłem testy porównawcze poniżej w celu rekomendacji numpy.corrcoef, głupio nie zdając sobie sprawy, że oryginalne pytanie już używa, corrcoefa w rzeczywistości pyta o pasowania wielomianowe wyższego rzędu. Dodałem rzeczywiste rozwiązanie do wielomianowego pytania r-kwadrat za pomocą modeli statystycznych i zostawiłem oryginalne testy porównawcze, które choć nie są na temat, są potencjalnie przydatne dla kogoś.

statsmodelsma możliwość bezpośredniego obliczenia r^2dopasowania wielomianu, oto 2 metody ...

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

Aby jeszcze bardziej skorzystać statsmodels, należy również przyjrzeć się dopasowanemu podsumowaniu modelu, które można wydrukować lub wyświetlić jako bogatą tabelę HTML w notatniku Jupyter / IPython. Obiekt wyników zapewnia dostęp do wielu przydatnych metryk statystycznych oprócz rsquared.

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

Poniżej znajduje się moja oryginalna odpowiedź, w której porównałem różne metody regresji liniowej r ^ 2 ...

Funkcja corrcoef użyta w pytaniu oblicza współczynnik korelacji r, tylko dla pojedynczej regresji liniowej, więc nie rozwiązuje problemu r^2pasowań wielomianowych wyższego rzędu. Jednak bez względu na to, co jest warte, doszedłem do wniosku, że w przypadku regresji liniowej jest to rzeczywiście najszybsza i najbardziej bezpośrednia metoda obliczania r.

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

Oto mój czas wynikający z porównania kilku metod dla 1000 losowych (x, y) punktów:

• Czysty Python (bezpośrednie robliczenia)
  • 1000 pętli, najlepiej 3: 1,59 ms na pętlę
• Numpy polyfit (dotyczy pasowań wielomianowych n-tego stopnia)
  • 1000 pętli, najlepiej 3: 326 µs na pętlę
• Numpy Manual (bezpośrednie robliczenia)
  • 10000 pętli, najlepiej 3: 62,1 µs na pętlę
• Numpy corrcoef (bezpośrednie robliczenia)
  • 10000 pętli, najlepiej 3: 56,6 µs na pętlę
• Scipy (regresja liniowa z rwyjściem)
  • 1000 pętli, najlepiej 3: 676 µs na pętlę
• Modele statystyk (mogą wykonywać wielomian n-tego stopnia i wiele innych dopasowań)
  • 1000 pętli, najlepiej 3: 422 µs na pętlę

Metoda corrcoef w niewielkim stopniu przewyższa obliczanie r ^ 2 „ręcznie” przy użyciu metod numpy. Jest> 5 razy szybszy niż metoda polyfit i ~ 12 razy szybszy niż scipy.linregress. Aby wzmocnić to, co robi dla ciebie numpy, jest 28 razy szybszy niż czysty Python. Nie jestem dobrze zorientowany w takich rzeczach, jak numba i pypy, więc ktoś inny musiałby wypełnić te luki, ale myślę, że jest to dla mnie bardzo przekonujące, że corrcoefjest to najlepsze narzędzie do obliczania rprostej regresji liniowej.

Oto mój kod do testów porównawczych. Skopiowałem i wkleiłem z notatnika Jupyter (trudno nie nazwać tego notatnikiem IPython ...), więc przepraszam, jeśli coś się zepsuło. Polecenie% timeit magic wymaga IPythona.

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-kwadrat to statystyka, która ma zastosowanie tylko do regresji liniowej.

Zasadniczo mierzy, jak dużą zmienność danych można wyjaśnić za pomocą regresji liniowej.

Zatem obliczasz „całkowitą sumę kwadratów”, która jest całkowitym kwadratem odchylenia każdej zmiennej wynikowej od ich średniej. . .

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

gdzie y_bar jest średnią z y.

Następnie obliczasz „regresyjną sumę kwadratów”, czyli o ile Twoje DOPASOWANE wartości różnią się od średniej

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

i znajdź stosunek tych dwóch.

Teraz wszystko, co musiałbyś zrobić, aby dopasować wielomian, to podłączyć y_hat pochodzi z tego modelu, ale nie jest dokładne nazywanie tego r-kwadrat.

Oto link, który znalazłem, który trochę do niego przemawia.

Artykuł Wikipedii na temat r-kwadratów sugeruje, że może być używany do ogólnego dopasowania modelu, a nie tylko do regresji liniowej.

Oto funkcja obliczająca ważone r-kwadrat za pomocą Pythona i Numpy (większość kodu pochodzi ze 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

Przykład:

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

wyjścia:

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

Odpowiada to formule ( lustro ):

gdzie f_i jest przewidywaną wartością z dopasowania, y_ {av} jest średnią z obserwowanych danych, y_i jest obserwowaną wartością danych. w_i to waga stosowana do każdego punktu danych, zwykle w_i = 1. SSE to suma kwadratów z powodu błędu, a SST to całkowita suma kwadratów.

Jeśli jesteś zainteresowany, kod w R: https://gist.github.com/dhimmel/588d64a73fa4fef02c8f ( mirror )

Oto bardzo prosta funkcja Pythona do obliczenia R ^ 2 na podstawie rzeczywistych i przewidywanych wartości, przy założeniu, że y i y_ to serie pand:

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)

Ze źródła scipy.stats.linregress. Używają metody średniej sumy kwadratów.

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

Możesz wykonać ten kod bezpośrednio, to znajdzie wielomian, a znajdziesz wartość R, którą możesz umieścić poniżej, jeśli potrzebujesz więcej wyjaśnień.

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)