Regresi Aljabar Kuadrat Terkecil Langkah-demi-Langkah Komputasi Aljabar Linier

22

Sebagai pendahuluan untuk pertanyaan tentang model linear-campuran dalam R, dan untuk berbagi sebagai referensi bagi para pecinta statistik pemula / menengah, saya memutuskan untuk memposting sebagai "gaya Tanya Jawab" yang independen, langkah-langkah yang terlibat dalam perhitungan "manual" dari koefisien dan nilai prediksi regresi linier sederhana.

Contohnya adalah dengan dataset R-built-in mtcars,, dan akan ditetapkan sebagai mil per galon yang dikonsumsi oleh kendaraan yang bertindak sebagai variabel independen, mundur dari bobot mobil (variabel kontinu), dan jumlah silinder sebagai faktor dengan tiga level (4, 6 atau 8) tanpa interaksi.

EDIT: Jika Anda tertarik dengan pertanyaan ini, Anda pasti akan menemukan jawaban yang terperinci dan memuaskan dalam posting ini oleh Matthew Drury di luar CV .

Antoni Parellada
sumber
Ketika Anda mengatakan "perhitungan manual", apa yang Anda cari? Ini relatif mudah untuk menunjukkan serangkaian langkah-langkah yang relatif sederhana untuk mendapatkan estimasi parameter dan sebagainya (melalui ortogonisasi Gram-Schmidt, misalnya, atau oleh operator SWEEP), tetapi bukan bagaimana R melakukan perhitungan secara internal; itu (dan sebagian besar paket statistik lainnya) menggunakan dekomposisi QR (dibahas dalam sejumlah posting di situs - pencarian pada dekomposisi QR menghasilkan sejumlah posting, beberapa di antaranya Anda mungkin mendapatkan nilai langsung dari)
Glen_b -Reinstate Monica
Iya nih. Saya percaya bahwa ini adalah alamat yang sangat baik dalam jawaban oleh MD. Saya mungkin harus mengedit posting saya, mungkin menekankan pendekatan geometris di belakang jawaban saya - ruang kolom, matriks proyeksi ...
Antoni Parellada
Ya! @Matthew Drury Apakah Anda ingin saya menghapus baris itu di OP, atau memperbarui tautan?
Antoni Parellada
1
Tidak yakin apakah Anda memiliki tautan ini, tetapi ini terkait erat, dan saya sangat suka jawaban JM. stats.stackexchange.com/questions/1829/…
Haitao Du

Jawaban:

51

Catatan : Saya telah memposting versi yang diperluas dari jawaban ini di situs web saya .

Apakah Anda akan mempertimbangkan untuk mengirim jawaban yang mirip dengan mesin R yang sebenarnya?

Yakin! Kita pergi ke lubang kelinci.

Lapisan pertama adalah lm, antarmuka terkena programmer R. Anda dapat melihat sumbernya hanya dengan mengetik lmdi konsol R. Mayoritas (seperti sebagian besar kode tingkat produksi) sibuk memeriksa input, mengatur atribut objek, dan melempar kesalahan; tapi garis ini menonjol

lm.fit(x, y, offset = offset, singular.ok = singular.ok, 
                ...)

lm.fitadalah fungsi R lain, Anda bisa menyebutnya sendiri. Sementara dengan lmmudah bekerja dengan formula dan bingkai data, membutuhkan lm.fitmatriks, jadi itu satu tingkat abstraksi dihapus. Memeriksa sumber untuk lm.fit, lebih banyak pekerjaan, dan baris yang sangat menarik berikut

z <- .Call(C_Cdqrls, x, y, tol, FALSE)

Sekarang kita menuju suatu tempat. .Calladalah cara R memanggil ke kode C. Ada fungsi C, C_Cdqrls di sumber R di suatu tempat, dan kita perlu menemukannya. Ini dia .

Melihat fungsi C, sekali lagi, kami menemukan sebagian besar batas memeriksa, pembersihan kesalahan, dan pekerjaan yang sibuk. Tetapi garis ini berbeda

F77_CALL(dqrls)(REAL(qr), &n, &p, REAL(y), &ny, &rtol,
        REAL(coefficients), REAL(residuals), REAL(effects),
        &rank, INTEGER(pivot), REAL(qraux), work);

Jadi sekarang kita menggunakan bahasa ketiga kita, R telah memanggil C yang memanggil fortran. Ini kode fortran .

Komentar pertama menceritakan semuanya

c     dqrfit is a subroutine to compute least squares solutions
c     to the system
c
c     (1)               x * b = y

(Menariknya, sepertinya nama rutin ini diubah di beberapa titik, tetapi seseorang lupa memperbarui komentar). Jadi kita akhirnya pada titik di mana kita dapat melakukan beberapa aljabar linier, dan benar-benar menyelesaikan sistem persamaan. Ini adalah hal yang benar-benar bagus untuk fortran, yang menjelaskan mengapa kami melewati begitu banyak lapisan untuk sampai ke sini.

Komentar tersebut juga menjelaskan apa yang akan dilakukan kode

c     on return
c
c        x      contains the output array from dqrdc2.
c               namely the qr decomposition of x stored in
c               compact form.

Jadi fortran akan memecahkan sistem dengan mencari dekomposisi.QR

Hal pertama yang terjadi, dan yang paling penting, adalah

call dqrdc2(x,n,n,p,tol,k,qraux,jpvt,work)

Ini memanggil fungsi fortran dqrdc2pada matriks input kami x. Apa ini?

 c     dqrfit uses the linpack routines dqrdc and dqrsl.

Jadi kami akhirnya berhasil linpack . Linpack adalah perpustakaan aljabar linear fortran yang telah ada sejak tahun 70-an. Aljabar linear yang paling serius akhirnya menemukan cara untuk linpack. Dalam kasus kami, kami menggunakan fungsi dqrdc2

c     dqrdc2 uses householder transformations to compute the qr
c     factorization of an n by p matrix x.

Di sinilah pekerjaan yang sebenarnya dilakukan. Butuh satu hari penuh yang baik bagi saya untuk mencari tahu apa yang dilakukan kode ini, levelnya serendah itu. Tetapi secara umum, kami memiliki matriks X dan kami ingin memfaktorkannya ke dalam produk mana Q adalah matriks ortogonal dan R adalah matriks segitiga atas. Ini adalah hal yang cerdas untuk dilakukan, karena begitu Anda memiliki Q dan R, Anda dapat menyelesaikan persamaan linear untuk regresiX=QRQRQR

XtXβ=XtY

sangat mudah. Memang

XtX=RtQtQR=RtR

jadi seluruh sistem menjadi

RtRβ=RtQty

tapi adalah segitiga atas dan memiliki peringkat yang sama seperti X t X , sehingga selama masalah kita baik yang ditimbulkan, itu adalah peringkat penuh, dan kami mungkin juga hanya memecahkan sistem berkurangRXtX

Rβ=Qty

Tapi ini hal yang luar biasa. adalah segitiga atas, jadi persamaan linear terakhir di sini adalah adil , jadi penyelesaian untuk β n adalah sepele. Anda kemudian bisa naik baris, satu per satu, dan gantikan di βRconstant * beta_n = constantβnβ s yang sudah Anda ketahui, setiap kali mendapatkan persamaan linear satu variabel sederhana untuk dipecahkan. Jadi, setelah Anda memiliki dan R , semuanya runtuh menjadi apa yang disebut substitusi mundur , yang mudah. Anda dapat membaca tentang ini secara lebih rinci di sini , di mana contoh kecil eksplisit sepenuhnya dikerjakan.QR

Matthew Drury
sumber
4
Ini adalah esai pendek matematika / coding paling menyenangkan yang bisa dibayangkan. Saya tidak tahu apa-apa tentang pengkodean, tetapi "tur" Anda melalui nyali fungsi R yang tampaknya tidak berbahaya benar-benar membuka mata. Tulisan yang sangat baik! Sejak "baik hati" melakukan trik ... Bisakah Anda berbaik hati mempertimbangkan satu ini sebagai tantangan terkait? :-)
Antoni Parellada
6
+1 Saya belum pernah melihat ini sebelumnya, ringkasan yang bagus. Hanya untuk menambahkan sedikit informasi jika @Antoni tidak terbiasa dengan transformasi Householder; itu pada dasarnya adalah transformasi linier yang memungkinkan Anda untuk menghilangkan satu bagian dari matriks R yang Anda coba capai tanpa memperbaiki bagian yang sudah Anda tangani (selama Anda melakukannya dengan urutan yang benar), menjadikannya ideal untuk mentransformasikan matriks ke bentuk segitiga atas (Rotasi Givens melakukan pekerjaan serupa dan mungkin lebih mudah divisualisasikan, tetapi sedikit lebih lambat). Saat Anda membangun R, Anda harus sekaligus membangun Q
Glen_b -Reinstate Monica
2
Matthew (+1), saya sarankan Anda memulai atau mengakhiri posting Anda dengan tautan ke madrury write-up Anda yang jauh lebih rinci .
Amuba kata Reinstate Monica
3
-1 untuk keluar dan tidak turun ke kode mesin.
S. Kolassa - Pasang kembali Monica
3
(Maaf, hanya bercanda ;-)
S. Kolassa - Reinstate Monica
8

Perhitungan langkah demi langkah aktual dalam R dijelaskan dengan indah dalam jawaban oleh Matthew Drury di utas yang sama ini. Dalam jawaban ini saya ingin berjalan melalui proses membuktikan kepada diri sendiri bahwa hasil dalam R dengan contoh sederhana dapat dicapai mengikuti aljabar linear proyeksi ke ruang kolom, dan konsep kesalahan tegak lurus (titik produk), diilustrasikan dalam posting yang berbeda , dan dijelaskan dengan baik oleh Dr. Strang di Aljabar Linier dan Penerapannya , dan mudah diakses di sini .

β

mhalg=sayantercehalt(cyl=4)+β1wesayaght+D1sayantercehalt(cyl=6)+D2sayantercehalt(cyl=8)[]

D1D2X

attach(mtcars)    
x1 <- wt

    x2 <- cyl; x2[x2==4] <- 1; x2[!x2==1] <-0

    x3 <- cyl; x3[x3==6] <- 1; x3[!x3==1] <-0

    x4 <- cyl; x4[x4==8] <- 1; x4[!x4==1] <-0

    X <- cbind(x1, x2, x3, x4)
    colnames(X) <-c('wt','4cyl', '6cyl', '8cyl')

head(X)
        wt 4cyl 6cyl 8cyl
[1,] 2.620    0    1    0
[2,] 2.875    0    1    0
[3,] 2.320    1    0    0
[4,] 3.215    0    1    0
[5,] 3.440    0    0    1
[6,] 3.460    0    1    0

[]lm

βPrHaijM.Sebuahtrsayax=(XTX)-1XT[PrHaijM.Sebuahtrsayax][y]=[RegrCHaiefs](XTX)-1XTy=β

X_tr_X_inv <- solve(t(X) %*% X)    
Proj_M <- X_tr_X_inv %*% t(X)
Proj_M %*% mpg

          [,1]
wt   -3.205613
4cyl 33.990794
6cyl 29.735212
8cyl 27.919934

Identik dengan: coef(lm(mpg ~ wt + as.factor(cyl)-1)).

HSebuahtM.Sebuahtrsayax=X(XTX)-1XT

HAT <- X %*% X_tr_X_inv %*% t(X)

y^X(XTX)-1XTy , dalam hal ini:, y_hat <- HAT %*% mpgyang memberikan nilai identik ke:

cyl <- as.factor(cyl); OLS <- lm(mpg ~ wt + cyl); predict(OLS):

y_hat <- as.numeric(y_hat)
predicted <- as.numeric(predict(OLS))
all.equal(y_hat,predicted)
[1] TRUE
Antoni Parellada
sumber
1
Secara umum, dalam komputasi numerik, saya percaya yang terbaik adalah menyelesaikan persamaan linear daripada menghitung matriks invers. Jadi, saya pikir beta = solve(t(X) %*% X, t(X) %*% y)dalam praktiknya lebih akurat daripada solve(t(X) %*% X) %*% t(X) %*% y.
Matthew Drury
R tidak melakukannya seperti itu - ia menggunakan dekomposisi QR. Jika Anda akan menggambarkan algoritma yang digunakan, pada komputer saya ragu ada yang menggunakan yang Anda tunjukkan.
Reinstate Monica - G. Simpson
Tidak setelah algoritma, hanya mencoba memahami dasar-dasar aljabar linier.
Antoni Parellada
@AntoniParellada Bahkan dalam kasus itu, saya masih menemukan pemikiran dalam hal persamaan linear yang lebih bersinar dalam banyak situasi.
Matthew Drury
1
Mengingat hubungan perifer dari utas ini dengan tujuan situs kami, sambil melihat nilai dalam menggambarkan penggunaan Runtuk perhitungan penting, saya ingin menyarankan agar Anda mempertimbangkan untuk mengubahnya menjadi kontribusi untuk blog kami.
whuber