ponderada - ¿Cómo puedo calcular r-squared usando Python y Numpy?

graficas estadisticas python

Question

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 que mejor se ajusta a Excel, y el valor r-cuadrado que calcula. Usando esto, sé que estoy calculando r-cuadrado correctamente para el mejor ajuste lineal (grado igual a 1). Sin embargo, mi función no funciona para polinomios con un grado mayor que 1.

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

Aquí está 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

Answer 1

Aquí hay una función para calcular el ponderado r-cuadrado 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

productos:

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

Esto corresponde a la formula ( mirror ):

con f_i es el valor previsto 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 a error y SST es la suma total de cuadrados.

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

Answer 2

De la documentación de numpy.polyfit , es apropiado 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 necesitas calcular el R-cuadrado para ese ajuste. La página de wikipedia sobre regresión lineal brinda detalles completos. Usted está interesado en R ^ 2, que puede calcular de varias maneras, siendo la más fácil probablemente

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''s, 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 manera más ordenada de calcular tu R-squared, 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

Answer 3

Desde yanl (yet-another-library) sklearn.metrics tiene una función r2_square ;

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

Answer 4

El artículo de wikipedia sobre en.wikipedia.org/wiki/Coefficient_of_determination sugiere que se puede usar para el ajuste general del modelo en lugar de solo la regresión lineal.

Answer 5

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

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 6

Originalmente publiqué los puntos de referencia a continuación con el propósito de recomendar numpy.corrcoef , neciamente sin darme cuenta de que la pregunta original ya usa corrcoef y de hecho estaba preguntando acerca de los ajustes polinómicos de orden superior. He agregado una solución real a la pregunta polinomial r-squared usando modelos de estadísticas, y dejé los puntos de referencia originales, que aunque fuera del tema, son potencialmente útiles para alguien.

statsmodels tiene la capacidad de calcular el r^2 de un ajuste polinómico directamente, aquí hay 2 métodos ...

import statsmodels.api as sm
import stasmodels.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

Para aprovechar aún más los statsmodels de statsmodels , también se debe mirar el resumen del modelo ajustado, que se puede imprimir o mostrar como una tabla HTML rica 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 única, por lo que no aborda la cuestión de r^2 para ajustes de polinomios de orden superior. Sin embargo, por lo que vale, he llegado a encontrar que para la regresión lineal, de hecho es el método más rápido y directo para calcular r .

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

Estos fueron mis resultados de tiempo al comparar un grupo de métodos para 1000 puntos aleatorios (x, y):

• Pure Python (cálculo directo r )
  • 1000 loops, lo mejor de 3: 1.59 ms por ciclo
• Numpty polyfit (aplicable a los ajustes polinomiales de grado n-ésimo)
  • 1000 loops, lo mejor de 3: 326 μs por ciclo
• Numpy Manual (cálculo directo r )
  • 10000 bucles, lo mejor de 3: 62.1 μs por ciclo
• Numex corrcoef (cálculo directo r )
  • 10000 loops, lo mejor de 3: 56.6 μs por ciclo
• Scipy (regresión lineal con r como salida)
  • 1000 bucles, lo mejor de 3: 676 μs por ciclo
• Statsmodels (puede hacer un polinomio de grado n y muchos otros ajustes)
  • 1000 bucles, lo mejor de 3: 422 μs por ciclo

El método de corrcoef supera ampliamente el cálculo de r ^ 2 "manualmente" usando métodos numpy. 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 soy 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 la corrcoef es la mejor herramienta para calcular r para una regresión lineal simple.

Aquí está mi código de evaluación comparativa. Copié y pegué de un cuaderno de Jupyter (difícil de no llamarlo un cuaderno de IPython ...), 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)
    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)

Answer 7

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

Básicamente, mide la cantidad de variación en sus datos que se puede explicar mediante la regresión lineal.

Entonces, se calcula la "Suma total de cuadrados", que es la desviación cuadrática total de cada una de las variables de resultado de su media. . .

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

donde y_bar es la media de las y.

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

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

y encuentra la proporción de esos dos.

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

Here hay un enlace que encontré que habla un poco.

Answer 8

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

scipy.stats.stats.linregress

es decir

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

como en la respuesta de @Adam Marples.