Bagaimana cara menghitung r-squared menggunakan Python dan Numpy?

Saya menggunakan Python dan Numpy untuk menghitung polinomial yang paling sesuai dari derajat arbitrer. Saya memberikan daftar nilai x, nilai y, dan derajat polinomial yang ingin saya paskan (linier, kuadrat, dll.).

Ini banyak berfungsi, tetapi saya juga ingin menghitung r (koefisien korelasi) dan r-kuadrat (koefisien determinasi). Saya membandingkan hasil saya dengan kapabilitas garis tren paling sesuai Excel, dan nilai r-squared yang dihitungnya. Dengan menggunakan ini, saya tahu saya menghitung r-kuadrat dengan benar untuk linier terbaik (derajat sama dengan 1). Namun, fungsi saya tidak berfungsi untuk polinomial dengan derajat lebih dari 1.

Excel mampu melakukan ini. Bagaimana cara menghitung r-kuadrat untuk polinomial orde tinggi menggunakan Numpy?

Inilah fungsi saya:

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

Dari dokumentasi numpy.polyfit , itu adalah regresi linier yang pas. Secara khusus, numpy.polyfit dengan derajat 'd' cocok dengan regresi linier dengan fungsi mean

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

Jadi, Anda hanya perlu menghitung R-kuadrat untuk kecocokan itu. Halaman wikipedia tentang regresi linier memberikan detail lengkap. Anda tertarik pada R ^ 2 yang dapat Anda hitung dengan beberapa cara, mungkin yang paling mudah

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

Di mana saya menggunakan 'y_bar' untuk mean dari y, dan 'y_ihat' untuk menjadi nilai yang cocok untuk setiap poin.

Saya tidak terlalu paham dengan numpy (saya biasanya bekerja di R), jadi mungkin ada cara yang lebih rapi untuk menghitung R-kuadrat Anda, tetapi yang berikut ini seharusnya benar

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

Balasan yang sangat terlambat, tetapi kalau-kalau seseorang membutuhkan fungsi siap untuk ini:

scipy.stats.linregress

yaitu

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

seperti dalam jawaban @Adam Marples.

Dari yanl (yet-another-library) sklearn.metricsmemiliki r2_scorefungsi;

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

Saya telah berhasil menggunakan ini, di mana x dan y mirip dengan array.

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

Saya awalnya memposting tolok ukur di bawah ini dengan tujuan merekomendasikan numpy.corrcoef, dengan bodohnya tidak menyadari bahwa pertanyaan asli sudah digunakan corrcoefdan sebenarnya menanyakan tentang kecocokan polinomial tingkat tinggi. Saya telah menambahkan solusi aktual untuk pertanyaan polinomial r-squared menggunakan statsmodels, dan saya telah meninggalkan tolok ukur asli, yang sementara di luar topik, berpotensi berguna bagi seseorang.

statsmodelsmemiliki kemampuan untuk menghitung r^2kesesuaian polinom secara langsung, berikut adalah 2 metode ...

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

Untuk memanfaatkan lebih jauh statsmodels, kita juga harus melihat ringkasan model yang dipasang, yang dapat dicetak atau ditampilkan sebagai tabel HTML kaya di notebook Jupyter / IPython. Objek hasil menyediakan akses ke banyak metrik statistik yang berguna sebagai tambahan rsquared.

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

Di bawah ini adalah Jawaban asli saya di mana saya membandingkan berbagai metode regresi linier r ^ 2 ...

Fungsi koreksi yang digunakan dalam Pertanyaan menghitung koefisien korelasi r, hanya untuk satu regresi linier, sehingga tidak menjawab pertanyaan tentang r^2kesesuaian polinomial orde tinggi. Namun, untuk apa nilainya, saya telah menemukan bahwa untuk regresi linier, itu memang metode penghitungan tercepat dan paling langsung r.

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

Ini adalah hasil waktu saya dari membandingkan sekumpulan metode untuk 1000 poin acak (x, y):

• Python murni ( rperhitungan langsung )
  • 1000 loop, terbaik 3: 1,59 ms per loop
• Poliit numpy (berlaku untuk kesesuaian polinomial derajat ke-n)
  • 1000 loop, terbaik 3: 326 µs per loop
• Numpy Manual ( rperhitungan langsung )
  • 10000 loop, terbaik 3: 62,1 µs per loop
• Numpy corrcoef ( rpenghitungan langsung )
  • 10000 loop, terbaik 3: 56,6 µs per loop
• Scipy (regresi linier dengan rsebagai keluaran)
  • 1000 loop, terbaik 3: 676 µs per loop
• Statsmodels (dapat melakukan polinomial derajat ke-n dan banyak kecocokan lainnya)
  • 1000 loop, terbaik 3: 422 µs per loop

Metode corrcoef mengalahkan penghitungan r ^ 2 "secara manual" menggunakan metode numpy. Ini> 5X lebih cepat dari metode polyfit dan ~ 12X lebih cepat dari scipy.linregress. Hanya untuk memperkuat apa yang numpy lakukan untuk Anda, ini 28X lebih cepat dari python murni. Saya tidak berpengalaman dalam hal-hal seperti numba dan pypy, jadi orang lain harus mengisi celah itu, tetapi saya pikir ini cukup meyakinkan bagi saya bahwa corrcoefini adalah alat terbaik untuk menghitung rregresi linier sederhana.

Ini kode pembandingan saya. Saya menyalin-tempel dari Notebook Jupyter (sulit untuk tidak menyebutnya sebagai Notebook IPython ...), jadi saya minta maaf jika ada yang rusak di jalan. Perintah sihir% timeit membutuhkan 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-squared adalah statistik yang hanya berlaku untuk regresi linier.

Pada dasarnya, ini mengukur seberapa banyak variasi dalam data Anda dapat dijelaskan dengan regresi linier.

Jadi, Anda menghitung "Jumlah Total Kuadrat", yang merupakan total deviasi kuadrat dari setiap variabel hasil Anda dari rata-ratanya. . .

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

dimana y_bar adalah mean dari y.

Kemudian, Anda menghitung "jumlah regresi kuadrat", yaitu seberapa besar nilai FITTED Anda berbeda dari rata-rata

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

dan temukan rasio keduanya.

Sekarang, yang harus Anda lakukan untuk pemasangan polinomial adalah menyambungkan y_hat dari model itu, tetapi tidak akurat untuk menyebutnya r-squared.

Ini adalah tautan yang saya temukan yang menjelaskannya sedikit.

Artikel wikipedia tentang r-squareds menunjukkan bahwa ini dapat digunakan untuk pemasangan model umum daripada hanya regresi linier.

Berikut adalah fungsi untuk menghitung weighted r-squared dengan Python dan Numpy (sebagian besar kode berasal dari 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

Contoh:

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

keluaran:

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

Ini sesuai dengan rumus ( cermin ):

dengan f_i adalah nilai prediksi dari fit, y_ {av} adalah mean dari data yang diamati y_i adalah nilai data yang diamati. w_i adalah pembobotan yang diterapkan pada setiap titik data, biasanya w_i = 1. SSE adalah jumlah kuadrat karena kesalahan dan SST adalah jumlah total kuadrat.

Jika tertarik, kodenya di R: https://gist.github.com/dhimmel/588d64a73fa4fef02c8f ( mirror )

Berikut adalah fungsi python yang sangat sederhana untuk menghitung R ^ 2 dari nilai aktual dan prediksi dengan asumsi y dan y_hat adalah seri 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)

Dari sumber scipy.stats.linregress. Mereka menggunakan metode jumlah rata-rata kuadrat.

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

Anda dapat menjalankan kode ini secara langsung, ini akan menemukan Anda polinomialnya, dan akan menemukan nilai-R untuk Anda yang dapat Anda beri komentar di bawah jika Anda memerlukan penjelasan lebih lanjut.

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)