Pierwotnie zamieściłem testy porównawcze poniżej w celu rekomendacji numpy.corrcoef, głupio nie zdając sobie sprawy, że oryginalne pytanie już używa, corrcoefa w rzeczywistości pyta o pasowania wielomianowe wyższego rzędu. Dodałem rzeczywiste rozwiązanie do wielomianowego pytania r-kwadrat za pomocą modeli statystycznych i zostawiłem oryginalne testy porównawcze, które choć nie są na temat, są potencjalnie przydatne dla kogoś.
statsmodelsma możliwość bezpośredniego obliczenia r^2dopasowania wielomianu, oto 2 metody ...
import statsmodels.api as sm
import statsmodels.formula.api as smf
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
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
Aby jeszcze bardziej skorzystać statsmodels, należy również przyjrzeć się dopasowanemu podsumowaniu modelu, które można wydrukować lub wyświetlić jako bogatą tabelę HTML w notatniku Jupyter / IPython. Obiekt wyników zapewnia dostęp do wielu przydatnych metryk statystycznych oprócz rsquared.
model = sm.OLS(y, xpoly)
results = model.fit()
results.summary()
Poniżej znajduje się moja oryginalna odpowiedź, w której porównałem różne metody regresji liniowej r ^ 2 ...
Funkcja corrcoef użyta w pytaniu oblicza współczynnik korelacji r, tylko dla pojedynczej regresji liniowej, więc nie rozwiązuje problemu r^2pasowań wielomianowych wyższego rzędu. Jednak bez względu na to, co jest warte, doszedłem do wniosku, że w przypadku regresji liniowej jest to rzeczywiście najszybsza i najbardziej bezpośrednia metoda obliczania r.
def get_r2_numpy_corrcoef(x, y):
return np.corrcoef(x, y)[0, 1]**2
Oto mój czas wynikający z porównania kilku metod dla 1000 losowych (x, y) punktów:
- Czysty Python (bezpośrednie
robliczenia)
- 1000 pętli, najlepiej 3: 1,59 ms na pętlę
- Numpy polyfit (dotyczy pasowań wielomianowych n-tego stopnia)
- 1000 pętli, najlepiej 3: 326 µs na pętlę
- Numpy Manual (bezpośrednie
robliczenia)
- 10000 pętli, najlepiej 3: 62,1 µs na pętlę
- Numpy corrcoef (bezpośrednie
robliczenia)
- 10000 pętli, najlepiej 3: 56,6 µs na pętlę
- Scipy (regresja liniowa z
rwyjściem)
- 1000 pętli, najlepiej 3: 676 µs na pętlę
- Modele statystyk (mogą wykonywać wielomian n-tego stopnia i wiele innych dopasowań)
- 1000 pętli, najlepiej 3: 422 µs na pętlę
Metoda corrcoef w niewielkim stopniu przewyższa obliczanie r ^ 2 „ręcznie” przy użyciu metod numpy. Jest> 5 razy szybszy niż metoda polyfit i ~ 12 razy szybszy niż scipy.linregress. Aby wzmocnić to, co robi dla ciebie numpy, jest 28 razy szybszy niż czysty Python. Nie jestem dobrze zorientowany w takich rzeczach, jak numba i pypy, więc ktoś inny musiałby wypełnić te luki, ale myślę, że jest to dla mnie bardzo przekonujące, że corrcoefjest to najlepsze narzędzie do obliczania rprostej regresji liniowej.
Oto mój kod do testów porównawczych. Skopiowałem i wkleiłem z notatnika Jupyter (trudno nie nazwać tego notatnikiem IPython ...), więc przepraszam, jeśli coś się zepsuło. Polecenie% timeit magic wymaga IPythona.
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)
