Wie kann ich r-Quadrat mit Python und Numpy berechnen?

StackOverflow https://stackoverflow.com/questions/893657

Frage

Ich verwende Python und Numpy ein Best-Fit-Polynom beliebigen Grades zu berechnen. Ich gebe eine Liste der x-Werte, y-Werte und der Grad des Polynoms Ich möchte passen (linear, quadratisch, usw.).

So viel funktioniert, aber ich möchte auch r (Korrelationskoeffizient) und r-squared (Bestimmtheitsmaß) berechnen. Ich vergleiche meine Ergebnisse mit Excel Best-Fit-Trendlinie Fähigkeit und den r-Quadrat-Wert berechnet sie. Mit dieser, ich weiß, dass ich die Berechnung r-Quadrat richtig für linearen Best-Fit (Grad gleich 1). Allerdings ist meine Funktion nicht für Polynome mit Grad arbeiten größer als 1 ist.

Excel ist in der Lage, dies zu tun. Wie berechne ich r-Quadrat für Polynomen 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
War es hilfreich?

Lösung

Von der numpy.polyfit Dokumentation, es ist passend, lineare Regression. Insbesondere numpy.polyfit mit Grad ‚d‘ paßt eine lineare Regression mit der mittleren Funktion

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

So Sie nur das R-Quadrat für diesen Sitz berechnen müssen. Die Wikipedia-Seite auf linearen Regression vollständige Details gibt. Sie interessieren sich für R ^ 2, die Sie in ein paar Möglichkeiten berechnen kann, die easisest wahrscheinlich zu sein 

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

Wo ich ‚y_bar‘ für den Mittelwert des y der Verwendung und ‚y_ihat‘ für jeden Punkt der passende Wert sein.

Ich bin nicht sehr vertraut mit numpy (I in der Regel in R arbeiten), so gibt es wahrscheinlich einen aufgeräumter Weg, um Ihr R-Quadrat zu berechnen, aber die folgenden soll 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

Andere Tipps

Eine sehr späte Antwort, aber nur für den Fall, dass jemand braucht eine Funktion bereit, für diese:

scipy.stats.linregress

d. 

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

wie in @ Adam Marples Antwort.

Von yanl (noch-another-Bibliothek) sklearn.metrics eine r2_square Funktion; 

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

Ich habe dies erfolgreich verwendet wird, wobei x und y-Array-like. 

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

I posted ursprünglich das Benchmarks unten mit dem Ziel numpy.corrcoef zu empfehlen, zu realisieren, dummerweise nicht, dass die ursprüngliche Frage corrcoef bereits verwendet und war in der Tat über höheres Ordnung Polynomanpassungen zu fragen. Ich habe eine tatsächliche Lösung für das Polynom r-Quadrat-Frage mit statsmodels hinzugefügt, und ich habe die ursprüngliche Benchmarks links, die während Wegthema, sind potentiell nützlich für jemanden.

statsmodels die Fähigkeit hat, die r^2 eines Polynomfit direkt zu berechnen, hier sind zwei 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

Zur weiteren Vorteil statsmodels zu nehmen, sollte man auf das angepasste Modell Zusammenfassung sehen auch, die gedruckt werden können oder als Rich-HTML-Tabelle angezeigt in Jupyter / IPython Notebook. Die Ergebnisse Objekt ermöglicht den Zugriff auf viele nützliche statistische Kennzahlen zusätzlich zu rsquared. 

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

Im Folgenden finden Sie meine ursprüngliche Antwort, wo ich verschiedene lineare Regression r ^ 2 Methoden gebenchmarkt ...

Die corrcoef Funktion in der Frage verwendet berechnet der Korrelationskoeffizient, r, nur für eine einzige lineare Regression, so dass es nicht die Frage nach r^2 höheren Ordnung Polynomanpassungen nicht ansprechen. Aber für das, was es wert ist, ich bin gekommen, dass für die lineare Regression zu finden, es ist in der Tat die schnellste und direkteste Methode r der Berechnung. 

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

Das waren meine timeit Ergebnisse aus dem Vergleich eine Reihe von Methoden für 1000 random (x, y) Punkte:

  • Reiner Python (direkte r Berechnung)
    • 1000 Loops, am besten von 3: 1,59 ms pro Schleife
  • Numpy polyfit (für n-ten Grades Polynomanpassungen)
    • 1000 Loops, am besten von 3: 326 & mgr; s pro Schleife
  • Numpy Manual (direkte r Berechnung)
    • 10000 Schleifen, am besten von 3: 62,1 & mgr; s pro Schleife
  • Numpy corrcoef (direkte r Berechnung)
    • 10000 Schleifen, am besten von 3: 56,6 & mgr; s pro Schleife
  • Scipy (lineare Regression mit r als Ausgang)
    • 1000 Loops, am besten von 3: 676 & mgr; s pro Schleife
  • Statsmodels (kann tun n-ten Grades Polynom und viele andere Anfälle)
    • 1000 Loops, am besten von 3: 422 & mgr; s pro Schleife

Die corrcoef Verfahren schlagen knapp die r ^ 2 „manuell“ die Berechnung unter Verwendung von numpy Methoden. Es ist> 5-mal schneller als die polyfit Verfahren und ~ 12x schneller als die scipy.linregress. Nur um zu verstärken, was numpy für Dich tut, dann ist es 28X schneller als reiner Python. Ich bin nicht sehr versiert in Sachen wie numba und PyPy, so jemand anderes diese Lücken füllen müßte, aber ich denke, das viel ist überzeugend für mich, dass corrcoef ist das beste Werkzeug r für eine einfache lineare Regression für die Berechnung.

Hier ist mein Benchmarking-Code. Ich kopiere kleistert von einem Jupyter Notebook (schwer, nicht zu nennen es ein IPython Notebook ...), so dass ich entschuldige mich, wenn etwas auf dem Weg brach. Das% timeit Magie Befehl 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)

Wikipedia-Artikel über r-squareds legt nahe, dass es für die allgemeine Modell verwendet werden kann, passend eher als nur lineare Regression.

Hier ist eine Funktion zur Berechnung der gewichtet r-Quadrat mit Python und Numpy (die meisten der Code stammt aus 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 die Formel ( Spiegel ):

mit f_i ist der vorhergesagte Wert aus dem Sitz, y_ {av} ist der Mittelwert der beobachteten Daten y_i die beobachtete Datenwert. w_i ist die Gewichtung für jeden Datenpunkt angewendet wird, in der Regel w_i = 1. SSE ist die Summe der Quadrate aufgrund eines Fehlers und SST ist die Gesamtsumme der Quadrate.

Wenn Sie interessiert sind, wird der Code in R: https://gist.github.com/dhimmel/588d64a73fa4fef02c8f ( Spiegel )

R-Quadrat ist eine Statistik, die lineare Regression gilt nur.

Im Wesentlichen misst, wie viel Variation in Ihren Daten durch die lineare Regression erklärt werden.

So berechnen Sie die „Gesamtsumme der Quadrate“, die die Gesamtabweichung von jedem Ihrer Ergebnisvariablen von ihrem Mittelwert quadriert ist. . .

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

Dabei gilt y_bar ist der Mittelwert der y ist.

Dann berechnen Sie die „Regressionsquadratsumme“, die, wie viel Ihre FITTED Werte vom Mittelwert abweichen ist

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

und das Verhältnis der beiden letztgenannten finden.

Nun, alles, was Sie würde für eine Polynomfit tun muß, ist Plug in der y_hat des von diesem Modell, aber es ist nicht richtig, dass r-squared zu nennen.

Hier ist ein Link ich fand, dass es spricht ein wenig.

Lizenziert unter: CC-BY-SA mit Zuschreibung
Nicht verbunden mit StackOverflow
scroll top