originalmente publicado por debajo de los puntos de referencia con el fin de recomendar numpy.corrcoef, tontamente sin darse cuenta de que la pregunta original ya utiliza corrcoef y de hecho fue preguntar por ataques de orden superior polinomio. 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 para calcular la 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
Tomar ventaja adicional de statsmodels, también hay que mirar el modelo ajustado resumen, que se puede imprimir o mostrar 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 es mi respuesta original, en el que como punto de referencia diferentes R^2 métodos de regresión lineal ...
La función corrcoef utilizado en la pregunta calcula el coeficiente de correlación, r, sólo para un único regresión lineal, por lo que no aborda la cuestión de r^2 para ajustes polinómicos de orden superior. Sin embargo, por lo que vale la pena, he llegado a encontrar que para la regresión lineal, de hecho es el método más rápido y directo de calcular r.
def get_r2_numpy_corrcoef(x, y):
return np.corrcoef(x, y)[0, 1]**2
Estas fueron mis TimeIt resultados de la comparación de un montón de métodos para 1000 aleatorio (X, Y) puntos:
- puro Python (
r cálculo directo)
- 1000 bucles, mejor de 3: 1.59 ms por ciclo
- Numphy polyfit (aplicable a los ajustes polinomiales de grado n-0)
- 1000 bucles, mejor de 3: 326 mu s por bucle
- Numpy Manual (
r cálculo directo)
- 10.000 bucles, mejor de 3: 62,1 mu s por bucle
- Numex corrcoef (cálculo directo)
- 10000 bucles, mejor de 3: 56,6 mu s por bucle
- Scipy (regresión lineal con
r como una salida)
- 1000 bucles, mejor de 3: 676 mu s por bucle
- Statsmodels (pueden hacer n-ésimo grado del polinomio y muchos otros ataques)
- 1000 bucles, mejor de 3: 422 mu s por bucle
El método de corrcoef supera por mucho el cálculo de la r^2 "manualmente" utilizando los 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 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)
Nota: se utiliza el grado sólo en el cálculo de coeffs. –
tydok es correcto. Está calculando la correlación de x, y, y r-cuadrado para y = p_0 + p_1 * x. Consulte mi respuesta a continuación para obtener un código que debería funcionar. Si no te importa que pregunte, ¿cuál es tu objetivo final? ¿Estás haciendo la selección del modelo (elegir qué grado usar)? ¿O algo mas? – leif
@leif - La solicitud se reduce a "hacerlo como lo hace Excel". Me da la sensación de estas respuestas de que los usuarios pueden estar leyendo demasiado en el valor r-cuadrado cuando se utiliza una curva no lineal de mejor ajuste. Sin embargo, no soy un asistente de matemáticas, y esta es la funcionalidad solicitada. –