Comment calculer le r-carré en utilisant Python et Numpy?

J'utilise Python et Numpy pour calculer un polynôme de meilleur ajustement de degré arbitraire. Je passe une liste de valeurs x, de valeurs y et du degré du polynôme que je veux ajuster (linéaire, quadratique, etc.).

Cela fonctionne beaucoup, mais je veux aussi calculer r (coefficient de corrélation) et r-carré (coefficient de détermination). Je compare mes résultats avec la capacité de courbe de tendance optimale d'Excel et la valeur au carré qu'il calcule. En utilisant cela, je sais que je calcule correctement le r-carré pour le meilleur ajustement linéaire (degré égal à 1). Cependant, ma fonction ne fonctionne pas pour les polynômes avec un degré supérieur à 1.

Excel est capable de le faire. Comment calculer le r-carré pour les polynômes d'ordre supérieur à l'aide de Numpy?

Voici ma fonction:

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

D'après la documentation numpy.polyfit , il s'agit d'une régression linéaire. Plus précisément, numpy.polyfit avec le degré 'd' correspond à une régression linéaire avec la fonction moyenne

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

Il vous suffit donc de calculer le R-carré pour cet ajustement. La page wikipedia sur la régression linéaire donne tous les détails. Vous êtes intéressé par R ^ 2 que vous pouvez calculer de plusieurs manières, la plus simple étant probablement

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

Où j'utilise 'y_bar' pour la moyenne des y, et 'y_ihat' pour être la valeur d'ajustement pour chaque point.

Je ne suis pas très familier avec numpy (je travaille généralement en R), il y a donc probablement un moyen plus propre de calculer votre R-carré, mais ce qui suit devrait être correct

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

Une réponse très tardive, mais juste au cas où quelqu'un aurait besoin d'une fonction prête pour cela:

scipy.stats.linregress

c'est à dire

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

comme dans la réponse de @Adam Marples.

De yanl (encore une autre bibliothèque) sklearn.metricsa une r2_scorefonction;

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

Je l'ai utilisé avec succès, où x et y ressemblent à un tableau.

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

J'ai initialement publié les points de repère ci-dessous dans le but de recommander numpy.corrcoef, stupidement ne me rendant pas compte que la question d'origine utilise déjà corrcoefet posait en fait des questions sur les ajustements polynomiaux d'ordre supérieur. J'ai ajouté une solution réelle à la question polynomiale r-carré à l'aide de statsmodels, et j'ai laissé les repères d'origine, qui, bien que hors sujet, sont potentiellement utiles à quelqu'un.

statsmodelsa la capacité de calculer r^2directement l'ajustement polynomial, voici 2 méthodes ...

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

Pour en tirer davantage parti statsmodels, il convient également de consulter le résumé du modèle ajusté, qui peut être imprimé ou affiché sous forme de tableau HTML riche dans le notebook Jupyter / IPython. L'objet de résultats donne accès à de nombreuses mesures statistiques utiles en plus de rsquared.

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

Voici ma réponse originale où j'ai comparé diverses méthodes de régression linéaire r ^ 2 ...

La fonction corrcoef utilisée dans la Question calcule le coefficient de corrélation r, uniquement pour une seule régression linéaire, elle ne répond donc pas à la question des r^2ajustements polynomiaux d'ordre supérieur. Cependant, pour ce que cela vaut, j'en suis venu à trouver que pour la régression linéaire, c'est en effet la méthode de calcul la plus rapide et la plus directe r.

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

Ce sont mes résultats timeit en comparant un tas de méthodes pour 1000 points aléatoires (x, y):

• Python pur ( rcalcul direct )
  • 1000 boucles, meilleur de 3: 1,59 ms par boucle
• Numpy polyfit (applicable aux ajustements polynomiaux au n-ième degré)
  • 1000 boucles, meilleur de 3: 326 µs par boucle
• Manuel Numpy ( rcalcul direct )
  • 10000 boucles, meilleur de 3: 62,1 µs par boucle
• Numpy corrcoef ( rcalcul direct )
  • 10000 boucles, meilleur de 3: 56,6 µs par boucle
• Scipy (régression linéaire avec rcomme sortie)
  • 1000 boucles, meilleur de 3: 676 µs par boucle
• Statsmodels (peut faire des polynômes au n-ième degré et de nombreux autres ajustements)
  • 1000 boucles, meilleur de 3: 422 µs par boucle

La méthode corrcoef bat étroitement le calcul du r ^ 2 "manuellement" en utilisant les méthodes numpy. Il est> 5 fois plus rapide que la méthode polyfit et ~ 12 fois plus rapide que scipy.linregress. Juste pour renforcer ce que numpy fait pour vous, il est 28 fois plus rapide que le python pur. Je ne connais pas bien des choses comme numba et pypy, donc quelqu'un d'autre devrait combler ces lacunes, mais je pense que c'est très convaincant pour moi que corrcoefc'est le meilleur outil pour calculer rune simple régression linéaire.

Voici mon code d'analyse comparative. J'ai copié-collé à partir d'un notebook Jupyter (difficile de ne pas l'appeler un notebook IPython ...), alors je m'excuse si quelque chose s'est cassé en cours de route. La commande% timeit magic nécessite 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_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)

Le R au carré est une statistique qui ne s'applique qu'à la régression linéaire.

Essentiellement, il mesure la variation de vos données qui peut être expliquée par la régression linéaire.

Ainsi, vous calculez la «somme totale des carrés», qui est l'écart carré total de chacune de vos variables de résultat par rapport à leur moyenne. . .

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

où y_bar est la moyenne des y.

Ensuite, vous calculez la "somme des carrés de régression", qui correspond à la différence entre vos valeurs FITTED et la moyenne

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

et trouvez le rapport de ces deux.

Maintenant, tout ce que vous auriez à faire pour un ajustement polynomial est de brancher les y_hat de ce modèle, mais ce n'est pas exact d'appeler cela r-carré.

Voici un lien que j'ai trouvé qui en parle un peu.

L'article de wikipedia sur les r-squareds suggère qu'il peut être utilisé pour l'ajustement général du modèle plutôt que pour la régression linéaire.

Voici une fonction pour calculer le r-carré pondéré avec Python et Numpy (la plupart du code provient de 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

Exemple:

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

les sorties:

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

Cela correspond à la formule ( miroir ):

avec f_i est la valeur prédite de l'ajustement, y_ {av} est la moyenne des données observées y_i est la valeur des données observées. w_i est la pondération appliquée à chaque point de données, généralement w_i = 1. SSE est la somme des carrés due à l'erreur et SST est la somme totale des carrés.

Si vous êtes intéressé, le code en R: https://gist.github.com/dhimmel/588d64a73fa4fef02c8f ( miroir )

Voici une fonction python très simple pour calculer R ^ 2 à partir des valeurs réelles et prévues en supposant que y et y_qui sont des séries pandas:

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)

De la source scipy.stats.linregress. Ils utilisent la méthode de la somme moyenne des carrés.

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

Vous pouvez exécuter ce code directement, cela vous trouvera le polynôme, et vous trouvera la valeur R, vous pouvez mettre un commentaire ci-dessous si vous avez besoin de plus d'explications.

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)