¿Cómo calculo r cuadrado usando Python y Numpy?

Estoy usando Python y Numpy para calcular un polinomio de mejor ajuste de grado arbitrario. Paso una lista de valores x, valores y y el grado del polinomio que quiero ajustar (lineal, cuadrático, etc.).

Esto funciona, pero también quiero calcular r (coeficiente de correlación) y r-cuadrado (coeficiente de determinación). Estoy comparando mis resultados con la capacidad de línea de tendencia de mejor ajuste de Excel y el valor r cuadrado que calcula. Al usar esto, sé que estoy calculando r cuadrado correctamente para el mejor ajuste lineal (el grado es igual a 1). Sin embargo, mi función no funciona para polinomios con grado mayor que 1.

Excel puede hacer esto. ¿Cómo calculo r cuadrado para polinomios de orden superior usando Numpy?

Esta es mi función:

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

De la documentación de numpy.polyfit , se ajusta a la regresión lineal. Específicamente, numpy.polyfit con grado 'd' se ajusta a una regresión lineal con la función media

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

Entonces, solo necesita calcular el R cuadrado para ese ajuste. La página de wikipedia sobre regresión lineal ofrece todos los detalles. Está interesado en R ^ 2, que puede calcular de un par de formas, probablemente la más sencilla

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

Donde uso 'y_bar' para la media de las y, y 'y_ihat' para ser el valor de ajuste para cada punto.

No estoy muy familiarizado con numpy (normalmente trabajo en R), por lo que probablemente haya una forma más ordenada de calcular su R cuadrado, pero lo siguiente debería ser correcto

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

Una respuesta muy tardía, pero en caso de que alguien necesite una función lista para esto:

scipy.stats.linregress

es decir

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

como en la respuesta de @Adam Marples.

De yanl (otra biblioteca más) sklearn.metricstiene una r2_scorefunción;

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

He estado usando esto con éxito, donde xey son como una matriz.

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

Originalmente publiqué los puntos de referencia a continuación con el propósito de recomendar numpy.corrcoef, tontamente sin darme cuenta de que la pregunta original ya se usa corrcoefy, de hecho, preguntaba acerca de los ajustes polinomiales de orden superior. Agregué una solución real a la pregunta del polinomio r-cuadrado usando modelos de estadísticas, y dejé los puntos de referencia originales, que aunque están fuera del tema, son potencialmente útiles para alguien.

statsmodelstiene la capacidad de calcular el r^2ajuste de un polinomio directamente, aquí hay 2 métodos ...

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

Para aprovechar aún más statsmodels, también se debe mirar el resumen del modelo ajustado, que puede imprimirse o mostrarse como una tabla HTML enriquecida en el cuaderno Jupyter / IPython. El objeto de resultados proporciona acceso a muchas métricas estadísticas útiles además de rsquared.

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

A continuación se muestra mi respuesta original donde comparé varios métodos de regresión lineal r ^ 2 ...

La función corrcoef utilizada en la Pregunta calcula el coeficiente de correlación r, solo para una regresión lineal simple, por lo que no aborda la cuestión de r^2los ajustes de polinomios de orden superior. Sin embargo, por lo que vale, he llegado a encontrar que para la regresión lineal, es de hecho el método de cálculo más rápido y directo r.

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

Estos fueron los resultados de mi tiempo al comparar varios métodos para 1000 puntos aleatorios (x, y):

• Python puro ( rcálculo directo )
  • 1000 bucles, lo mejor de 3: 1,59 ms por bucle
• Numpy polyfit (aplicable a ajustes polinomiales de n-ésimo grado)
  • 1000 bucles, lo mejor de 3: 326 µs por bucle
• Numpy Manual ( rcálculo directo )
  • 10000 bucles, lo mejor de 3: 62,1 µs por bucle
• Numpy corrcoef ( rcálculo directo )
  • 10000 bucles, lo mejor de 3: 56,6 µs por bucle
• Scipy (regresión lineal con rcomo salida)
  • 1000 bucles, lo mejor de 3: 676 µs por bucle
• Statsmodels (puede hacer polinomios de n-ésimo grado y muchos otros ajustes)
  • 1000 bucles, lo mejor de 3: 422 µs por bucle

El método corrcoef supera por poco el cálculo de r ^ 2 "manualmente" utilizando métodos numerosos. Es> 5 veces más rápido que el método polyfit y ~ 12 veces más rápido que el scipy.linregress. Solo para reforzar lo que Numpy está haciendo por ti, es 28 veces más rápido que Python puro. No estoy muy versado en cosas como numba y pypy, por lo que alguien más tendría que llenar esos vacíos, pero creo que esto me convence bastante de que corrcoefes la mejor herramienta para calcular runa regresión lineal simple.

Aquí está mi código de evaluación comparativa. Copié y pegué de un Jupyter Notebook (es difícil no llamarlo un IPython Notebook ...), así que me disculpo si algo se rompió en el camino. El comando% timeit magic requiere 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)

R cuadrado es una estadística que solo se aplica a la regresión lineal.

Esencialmente, mide cuánta variación en sus datos se puede explicar mediante la regresión lineal.

Entonces, calcula la "Suma total de cuadrados", que es la desviación al cuadrado total de cada una de las variables de resultado de su media. . .

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

donde y_bar es la media de las y.

Luego, calcula la "suma de cuadrados de regresión", que es cuánto difieren sus valores FITTED de la media.

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

y encuentra la razón de esos dos.

Ahora, todo lo que tendría que hacer para un ajuste polinomial es conectar el y_que es de ese modelo, pero no es exacto llamarlo r cuadrado.

Aquí hay un enlace que encontré que habla un poco.

El artículo de wikipedia sobre r-cuadrados sugiere que se puede usar para el ajuste general del modelo en lugar de solo la regresión lineal.

Aquí hay una función para calcular el r-cuadrado ponderado con Python y Numpy (la mayor parte del código proviene 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

Ejemplo:

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

salidas:

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

Esto corresponde a la fórmula ( espejo ):

con f_i es el valor predicho del ajuste, y_ {av} es la media de los datos observados y_i es el valor de los datos observados. w_i es la ponderación aplicada a cada punto de datos, generalmente w_i = 1. SSE es la suma de cuadrados debido al error y SST es la suma total de cuadrados.

Si está interesado, el código en R: https://gist.github.com/dhimmel/588d64a73fa4fef02c8f ( espejo )

Aquí hay una función de Python muy simple para calcular R ^ 2 a partir de los valores reales y predichos, asumiendo que y e y_que son series de 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 fuente scipy.stats.linregress. Usan el método de suma de cuadrados promedio.

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

Puede ejecutar este código directamente, esto le encontrará el polinomio, y le encontrará el valor R, puede poner un comentario abajo si necesita más explicación.

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)