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

Question

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

Cela fonctionne, 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é d'Excel à calculer la meilleure ligne de tendance et la valeur r-carré qu'il calcule. Je sais ainsi que je calcule correctement le r-carré pour un meilleur ajustement linéaire (degré égal à 1). Cependant, ma fonction ne fonctionne pas pour les polynômes dont le degré est supérieur à 1.

Excel est en mesure 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

Answer 1

Une réponse très tardive, mais juste au cas où quelqu'un aurait besoin d'une fonction prête à l'emploi 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.

Answer 2

De yanl (yet-another-library) sklearn.metrics dispose d'un r2_score fonction ;

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

Answer 3

A partir de la numpy.polyfit la documentation, il s'agit d'une régression linéaire. Plus précisément, numpy.polyfit avec le degré 'd' ajuste une régression linéaire avec la fonction mean

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

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

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

J'utilise 'y_bar' pour la moyenne des y et 'y_ihat' pour la valeur d'ajustement de chaque point.

Je ne suis pas très familier avec numpy (je travaille habituellement avec R), il y a donc probablement une façon plus simple 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

Answer 4

J'ai utilisé cette méthode avec succès, où x et y sont des tableaux.

Note : pour la régression linéaire uniquement

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

Answer 5

À l'origine, j'ai publié les critères d'évaluation ci-dessous dans le but de recommander numpy.corrcoef , ne réalisant pas que la question originale utilise déjà le terme corrcoef et posait en fait une question sur les ajustements polynomiaux d'ordre supérieur. J'ai ajouté une solution réelle à la question du r-carré polynomial en utilisant statsmodels, et j'ai laissé les repères originaux qui, bien que hors sujet, sont potentiellement utiles à quelqu'un.

---

statsmodels a la capacité de calculer la r^2 d'un ajustement polynomial directement, 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 mieux tirer parti de la statsmodels Il convient également de consulter le résumé du modèle ajusté, qui peut être imprimé ou affiché sous la forme d'un tableau HTML riche dans le carnet Jupyter/IPython. L'objet results permet d'accéder à de nombreuses mesures statistiques utiles en plus de l'objet rsquared .

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

---

Vous trouverez ci-dessous ma réponse originale dans laquelle j'ai comparé différentes méthodes de régression linéaire r^2...

En corrcoef utilisée dans la Question calcule le coefficient de corrélation, r ne s'applique qu'à une seule régression linéaire, de sorte qu'elle ne répond pas à la question de savoir si l'on a besoin d'une régression linéaire ou d'une régression linéaire. r^2 pour les ajustements polynomiaux d'ordre supérieur. Cependant, pour ce que cela vaut, j'ai constaté que pour la régression linéaire, c'est en effet la méthode la plus rapide et la plus directe de calculer r .

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

Voici les résultats que j'ai obtenus en comparant un certain nombre de méthodes pour 1000 points aléatoires (x, y) :

• Pure Python (direct r calcul)
  • 1000 boucles, meilleur des 3 : 1,59 ms par boucle
• Numpy polyfit (applicable aux ajustements polynomiaux du troisième degré)
  • 1000 boucles, meilleur des 3 : 326 µs par boucle
• Manuel Numpy (direct) r calcul)
  • 10000 boucles, meilleur des 3 : 62,1 µs par boucle
• Numpy corrcoef (direct r calcul)
  • 10000 boucles, meilleur des 3 : 56,6 µs par boucle
• Scipy (régression linéaire avec r en sortie)
  • 1000 boucles, meilleur des 3 : 676 µs par boucle
• Statsmodels (peut effectuer des ajustements polynomiaux du troisième degré et de nombreux autres ajustements)
  • 1000 boucles, meilleur des 3 : 422 µs par boucle

La méthode corrcoef l'emporte de peu sur le calcul "manuel" de r^2 à l'aide des méthodes numpy. Elle est >5X plus rapide que la méthode polyfit et ~12X plus rapide que scipy.linregress. Juste pour renforcer ce que numpy fait pour vous, c'est 28X plus rapide que python pur. Je ne suis pas très au fait de choses comme numba et pypy, donc quelqu'un d'autre devra combler ces lacunes, mais je pense que c'est suffisamment convaincant pour moi que corrcoef est le meilleur outil pour calculer r pour une régression linéaire simple.

Voici mon code de référence. J'ai fait un copier-coller à partir d'un Notebook Jupyter (difficile de ne pas l'appeler Notebook IPython...), donc je m'excuse si quelque chose s'est cassé en chemin. La commande magique %timeit 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)

7/28/21 Résultats du benchmarking. (Python 3.7, numpy 1.19, scipy 1.6, statsmodels 0.12)

Python
2.41 ms ± 180 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Numpy polyfit
318 µs ± 44.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Numpy Manual
79.3 µs ± 4.05 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
Numpy corrcoef
83.8 µs ± 1.37 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
Scipy
221 µs ± 7.12 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Statsmodels
375 µs ± 3.63 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)