numpy.corrcoef私は元々 、元の質問がすでに使用されていることに気づかずcorrcoef、実際には高階多項式の適合について質問していたことを推奨する目的で、以下のベンチマークを投稿しました。statsmodelsを使用して、多項式のr-squared質問に実際のソリューションを追加しました。元のベンチマークは残しました。これは、トピックから外れていますが、誰かに役立つ可能性があります。
statsmodelsr^2多項式フィットのを直接計算する機能があります。ここに2つの方法があります。
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
をさらに活用するにはstatsmodels、Jupyter/IPythonノートブックでリッチHTMLテーブルとして印刷または表示できるフィットモデルの概要も確認する必要があります。結果オブジェクトは、に加えて多くの有用な統計メトリックへのアクセスを提供しますrsquared。
model = sm.OLS(y, xpoly)
results = model.fit()
results.summary()
以下は、さまざまな線形回帰r^2メソッドのベンチマークを行った元の回答です...
質問で使用されるcorrcoef関数は、単一の線形回帰に対してのみ相関係数を計算するため、高階多項式近似rの問題には対応していません。r^2ただし、その価値については、線形回帰の場合、実際に最も速く、最も直接的な計算方法であることがわかりましたr。
def get_r2_numpy_corrcoef(x, y):
return np.corrcoef(x, y)[0, 1]**2
これらは、1000個のランダム（x、y）ポイントの一連のメソッドを比較した結果です。
- 純粋なPython（直接
r計算）
- 1000ループ、ベスト3：ループあたり1.59ミリ秒
- Numpy polyfit（n次多項式近似に適用可能）
- 1000ループ、ベスト3：ループあたり326 µs
- Numpyマニュアル（直接
r計算）
- 10000ループ、ベスト3：ループあたり62.1 µs
- Numpy corrcoef（直接
r計算）
- 10000ループ、ベスト3：ループあたり56.6 µs
- Scipy（
r出力としての線形回帰）
- 1000ループ、ベスト3：ループあたり676 µs
- Statsmodels（n次多項式および他の多くの近似を実行できます）
- 1000ループ、ベスト3：ループあたり422 µs
corrcoefメソッドは、numpyメソッドを使用して「手動で」r^2を計算するよりもわずかに優れています。polyfitメソッドよりも5倍以上速く、scipy.linregressよりも12倍速くなります。numpyがあなたのために何をしているかを補強するために、それは純粋なpythonより28倍高速です。私はnumbaやpypyのようなものに精通していないので、他の誰かがそれらのギャップを埋める必要がありますが、これは単純な線形回帰corrcoefを計算するための最良のツールであると私には十分に説得力があると思います。r
これが私のベンチマークコードです。Jupyter Notebookからコピーして貼り付けたので（IPython Notebookとは呼ばないでください...）、途中で何かが壊れた場合はお詫びします。％timeitmagicコマンドには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ベンチマーク結果。（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)
