Wie berechne ich das R-Quadrat mit Python und Numpy?

Ich verwende Python und Numpy, um ein Best-Fit-Polynom beliebigen Grades zu berechnen. Ich übergebe eine Liste mit x-Werten, y-Werten und dem Grad des Polynoms, das ich anpassen möchte (linear, quadratisch usw.).

So viel funktioniert, aber ich möchte auch r (Korrelationskoeffizient) und r-Quadrat (Bestimmungskoeffizient) berechnen. Ich vergleiche meine Ergebnisse mit der am besten passenden Trendlinienfunktion von Excel und dem berechneten R-Quadrat-Wert. Auf diese Weise weiß ich, dass ich das r-Quadrat für die lineare Bestanpassung korrekt berechne (Grad gleich 1). Meine Funktion funktioniert jedoch nicht für Polynome mit einem Grad größer als 1.

Excel ist dazu in der Lage. Wie berechne ich das r-Quadrat für Polynome höherer Ordnung mit Numpy?

Hier ist meine Funktion:

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

In der Dokumentation zu numpy.polyfit wird die lineare Regression angepasst. Insbesondere passt numpy.polyfit mit dem Grad 'd' eine lineare Regression mit der mittleren Funktion an

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

Sie müssen also nur das R-Quadrat für diese Anpassung berechnen. Die Wikipedia-Seite zur linearen Regression enthält alle Details. Sie interessieren sich für R ^ 2, das Sie auf verschiedene Arten berechnen können, wobei es wahrscheinlich am einfachsten ist

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

Wobei ich 'y_bar' als Mittelwert der ys und 'y_ihat' als Anpassungswert für jeden Punkt verwende.

Ich bin mit Numpy nicht besonders vertraut (ich arbeite normalerweise in R), daher gibt es wahrscheinlich eine übersichtlichere Methode, um Ihr R-Quadrat zu berechnen, aber das Folgende sollte korrekt sein

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

Eine sehr späte Antwort, aber nur für den Fall, dass jemand eine fertige Funktion dafür benötigt:

scipy.stats.linregress

dh

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

wie in der Antwort von @Adam Marples.

From yanl (noch eine andere Bibliothek) sklearn.metricshat eine r2_scoreFunktion;

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

Ich habe dies erfolgreich verwendet, wobei x und y Array-ähnlich sind.

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

Ich habe ursprünglich die folgenden Benchmarks mit dem Ziel veröffentlicht, zu empfehlen numpy.corrcoef, dummerweise nicht zu erkennen, dass die ursprüngliche Frage bereits verwendet wird corrcoefund tatsächlich nach Polynomanpassungen höherer Ordnung gefragt hat. Ich habe eine tatsächliche Lösung für die Polynom-R-Quadrat-Frage mithilfe von Statistikmodellen hinzugefügt und die ursprünglichen Benchmarks beibehalten, die zwar nicht zum Thema gehören, aber möglicherweise für jemanden nützlich sind.

statsmodelshat die Fähigkeit, die r^2Polynomanpassung direkt zu berechnen , hier sind 2 Methoden ...

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

Um dies weiter zu nutzen statsmodels, sollte man sich auch die angepasste Modellzusammenfassung ansehen, die in Jupyter / IPython-Notizbuch gedruckt oder als umfangreiche HTML-Tabelle angezeigt werden kann. Das Ergebnisobjekt bietet zusätzlich Zugriff auf viele nützliche statistische Metriken rsquared.

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

Unten ist meine ursprüngliche Antwort, in der ich verschiedene lineare Regressionsmethoden verglichen habe ...

Die in der Frage verwendete Corrcoef- Funktion berechnet den Korrelationskoeffizienten rnur für eine einzelne lineare Regression, sodass die Frage r^2nach Polynomanpassungen höherer Ordnung nicht behandelt wird . Für das, was es wert ist, habe ich jedoch festgestellt, dass es für die lineare Regression tatsächlich die schnellste und direkteste Berechnungsmethode ist r.

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

Dies waren meine zeitlichen Ergebnisse aus dem Vergleich einer Reihe von Methoden für 1000 zufällige (x, y) Punkte:

• Pure Python (direkte rBerechnung)
  • 1000 Schleifen, am besten 3: 1,59 ms pro Schleife
• Numpy Polyfit (anwendbar auf Polynomanpassungen n-ten Grades)
  • 1000 Schleifen, am besten 3: 326 µs pro Schleife
• Numpy Manual (direkte rBerechnung)
  • 10000 Schleifen, am besten 3: 62,1 µs pro Schleife
• Numpy Corrcoef (direkte rBerechnung)
  • 10000 Schleifen, am besten 3: 56,6 µs pro Schleife
• Scipy (lineare Regression mit rals Ausgabe)
  • 1000 Schleifen, am besten 3: 676 µs pro Schleife
• Statistikmodelle (können Polynome n-ten Grades und viele andere Anpassungen ausführen)
  • 1000 Schleifen, am besten 3: 422 µs pro Schleife

Die Corrcoef-Methode schlägt die Berechnung des r ^ 2 "manuell" mit Numpy-Methoden knapp. Es ist> 5X schneller als die Polyfit-Methode und ~ 12X schneller als die scipy.linregress. Nur um zu verstärken, was Numpy für Sie tut, ist es 28-mal schneller als reines Python. Ich bin mit Dingen wie Numba und Pypy nicht vertraut, daher müsste jemand anderes diese Lücken füllen, aber ich denke, das überzeugt mich sehr, dass dies corrcoefdas beste Werkzeug für die Berechnung reiner einfachen linearen Regression ist.

Hier ist mein Benchmarking-Code. Ich habe von einem Jupyter-Notizbuch kopiert (es ist schwer, es nicht als IPython-Notizbuch zu bezeichnen ...), also entschuldige ich mich, wenn unterwegs etwas kaputt gegangen ist. Der Befehl% timeit magic erfordert 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)

Das R-Quadrat ist eine Statistik, die nur für die lineare Regression gilt.

Im Wesentlichen wird gemessen, wie stark die Variation Ihrer Daten durch die lineare Regression erklärt werden kann.

Sie berechnen also die "Gesamtsumme der Quadrate", die die gesamte quadratische Abweichung jeder Ihrer Ergebnisvariablen von ihrem Mittelwert darstellt. . .

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

Dabei ist y_bar der Mittelwert der ys.

Dann berechnen Sie die "Regressionssumme der Quadrate", dh wie stark sich Ihre FITTED-Werte vom Mittelwert unterscheiden

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

und finde das Verhältnis dieser beiden.

Alles, was Sie für eine Polynomanpassung tun müssten, ist das Einstecken des y_hat aus diesem Modell, aber es ist nicht genau, dieses r-Quadrat zu nennen.

Hier ist ein Link, den ich gefunden habe und der ein wenig dazu spricht.

Der Wikipedia-Artikel über R-Quadrate schlägt vor, dass er eher für die allgemeine Modellanpassung als nur für die lineare Regression verwendet werden kann.

Hier ist eine Funktion zum Berechnen des gewichteten R-Quadrats mit Python und Numpy (der größte Teil des Codes stammt von 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

Beispiel:

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

Ausgänge:

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

Dies entspricht der Formel ( Spiegel ):

mit f_i ist der vorhergesagte Wert aus der Anpassung, y_ {av} ist der Mittelwert der beobachteten Daten y_i ist der beobachtete Datenwert. w_i ist die Gewichtung, die auf jeden Datenpunkt angewendet wird, normalerweise w_i = 1. SSE ist die Summe der fehlerhaften Quadrate und SST ist die Gesamtsumme der Quadrate.

Bei Interesse den Code in R: https://gist.github.com/dhimmel/588d64a73fa4fef02c8f ( Spiegel )

Hier ist eine sehr einfache Python-Funktion zum Berechnen von R ^ 2 aus den tatsächlichen und vorhergesagten Werten unter der Annahme, dass y und y_pandas-Reihen sind:

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)

Aus der Quelle scipy.stats.linregress. Sie verwenden die Methode der durchschnittlichen Quadratsumme.

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

Sie können diesen Code direkt ausführen. Hier finden Sie das Polynom und den R-Wert. Wenn Sie weitere Erläuterungen benötigen, können Sie unten einen Kommentar abgeben.

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)