Perhitungan variabel kontras polinomial

11

Tolong beri saya ide bagaimana secara efisien mengode ulang variabel kategori (faktor) ke dalam set variabel kontras polinomial ortogonal.

Untuk banyak jenis variabel kontras (mis. Penyimpangan, simpel, Helmert, dll.) Pass-nya adalah:

  1. Tulis matriks koefisien kontras sesuai dengan jenisnya.
  2. Membalikkan atau menggeneralisasi-membaliknya untuk mendapatkan matriks kode.

Sebagai contoh:

Suppose there is 3-group factor and we want to recode it into a set of deviation  contrast variables.
The last group is treated as reference. Then the contrast coefficients matrix L is

         Group1 Group2 Group3
   var1   2/3   -1/3   -1/3
   var2  -1/3    2/3   -1/3

and ginv(L) is then the sought-for coding matrix

         var1 var2
  Group1   1    0
  Group2   0    1
  Group3  -1   -1

(We might also use inv(L) instead if we add a row for constant, equal to 1/3, at the head of L.)

Apakah ada cara yang sama atau serupa untuk mendapatkan variabel kontras polinomial? Jika ya, seperti apa matriks C dan bagaimana cara menyusunnya? Jika tidak, apa yang masih merupakan cara untuk menghitung variabel kontras polinom secara efisien (misalnya dengan aljabar matriks).

ttnphns
sumber
1
Saya melihat pertanyaan Anda setelah memverifikasi (kebetulan) bahwa hasil qr.qy()setuju dengan perhitungan manual qr.Q(qr(X))diikuti oleh Q%*%zpada posting saya. Saya benar-benar ingin tahu apakah saya dapat mengatakan sesuatu yang berbeda untuk menjawab pertanyaan Anda tanpa duplikasi. Saya benar-benar tidak ingin melakukan pekerjaan yang buruk ... Saya telah membaca cukup banyak posting Anda sehingga sangat menghormati Anda ... Jika saya menemukan cara untuk mengekspresikan konsep tanpa kode, secara konseptual hanya melalui aljabar linier, Saya akan kembali ke sana. Saya senang, bagaimanapun, bahwa Anda menemukan eksplorasi saya tentang masalah beberapa nilai. Salam hangat, Toni.
Antoni Parellada
@Antoni, terima kasih. Tujuan saya adalah untuk dapat mengkodekannya (dalam SPSS, dengan sintaksnya). Apakah mungkin untuk mengetahui cara kerja fungsi yang Anda sebutkan?
ttnphns

Jawaban:

5

Sebagai tambahan untuk posting saya sebelumnya pada topik ini saya ingin berbagi beberapa sementara (meskipun tidak lengkap) eksplorasi fungsi di belakang aljabar linier dan fungsi R terkait. Ini seharusnya menjadi pekerjaan yang sedang berjalan.

Bagian dari ketidakjelasan fungsi berkaitan dengan bentuk "kompak" dari dekomposisi Householder . Gagasan di balik dekomposisi Householder adalah untuk merefleksikan vektor melintasi sebuah hyperplane yang ditentukan oleh vektor-unit seperti pada diagram di bawah ini, tetapi memilih pesawat ini dengan cara yang bertujuan untuk memproyeksikan setiap vektor kolom dari matriks asli ke vektor satuan standar. Vektor norma-2 dinormalisasi dapat digunakan untuk menghitung transformasi Householder yang berbeda .u A e 1 1 u I - 2QRuAe11uI2uuTx

masukkan deskripsi gambar di sini

Proyeksi yang dihasilkan dapat dinyatakan sebagai

sign(xi=x1)×x[1000]+[x1x2x3xm]

Vektor mewakili perbedaan antara vektor kolom dalam matriks yang ingin kita dekomposisi dan vektor berhubungan dengan refleksi melintasi subruang atau "mirror" yang ditentukan oleh .vxAyu

Metode yang digunakan oleh LAPACK membebaskan kebutuhan untuk penyimpanan entri pertama di reflektor Householder dengan mengubahnya menjadi 's. Alih-alih menormalkan vektor ke dengan , itu hanya entri pertama yang dikonversi menjadi ; namun, vektor-vektor baru ini - sebut saja masih dapat digunakan sebagai vektor arah.1vuu=11w

Keindahan dari metode ini adalah bahwa mengingat bahwa dalam dekomposisi adalah segitiga atas, kita benar-benar dapat mengambil keuntungan dari elemen dalam di bawah diagonal untuk mengisi mereka dengan ini reflektor. Untungnya, entri-entri utama dalam vektor-vektor ini semuanya sama dengan , mencegah masalah dalam diagonal "yang disengketakan" dari matriks: mengetahui bahwa semuanya adalah mereka tidak perlu dimasukkan, dan dapat menghasilkan diagonal ke entri .RQR0Rw11R

Matriks "compact QR" dalam fungsi qr()$qrdapat dipahami secara kasar sebagai penambahan matriks dan matriks "storage" segitiga yang lebih rendah untuk reflektor yang "dimodifikasi".R

Proyeksi Householder masih akan memiliki bentuk , tetapi kami tidak akan bekerja dengan ( ), melainkan dengan vektor , yang hanya entri pertama yang dijamin , danI2uuTxux=1w1

(1)I2uuTx=I2wwwTwx=I2wwTw2x .

Orang akan berasumsi bahwa itu akan baik-baik saja untuk menyimpan ini reflektor di bawah diagonal atau tidak termasuk entri pertama dari , dan menyebutnya sehari. Namun, segala sesuatunya tidak pernah semudah ini. Alih-alih apa yang disimpan di bawah diagonal dalam adalah kombinasi dan koefisien dalam transformasi Householder dinyatakan sebagai (1), sehingga, mendefinisikan sebagai:wR1qr()$qrwtau

τ=wTw2=w2 , reflektor dapat diekspresikan sebagai . Vektor "reflektor" ini adalah yang disimpan tepat di bawah dalam apa yang disebut "compact ".reflectors=w/τRQR

Sekarang kita satu derajat menjauh dari vektor, dan entri pertama tidak lagi , Oleh karena itu output dari kehendak perlu menyertakan kunci untuk mengembalikan mereka karena kita bersikeras tidak termasuk masuknya pertama dari "reflektor" vektor untuk cocok semuanya . Jadi, apakah kita melihat nilai dalam output? Yah, tidak ada yang bisa diprediksi. Alih-alih dalam output (tempat kunci ini disimpan) kita menemukan .w1qr()qr()$qrτqr()$qrauxρ=reflectors22=wTwτ2/2

Jadi dibingkai dengan warna merah di bawah ini, kita melihat "reflektor" ( ), tidak termasuk entri pertama mereka.w/τ

masukkan deskripsi gambar di sini

Semua kode ada di sini , tetapi karena jawaban ini adalah tentang persimpangan kode dan aljabar linier, saya akan menempelkan hasilnya untuk memudahkan:


options(scipen=999)
set.seed(13)
(X = matrix(c(rnorm(16)), nrow=4, byrow=F))
           [,1]      [,2]       [,3]       [,4]
[1,]  0.5543269 1.1425261 -0.3653828 -1.3609845
[2,] -0.2802719 0.4155261  1.1051443 -1.8560272
[3,]  1.7751634 1.2295066 -1.0935940 -0.4398554
[4,]  0.1873201 0.2366797  0.4618709 -0.1939469

Sekarang saya menulis fungsinya House()sebagai berikut:

   House = function(A){
    Q = diag(nrow(A))
    reflectors = matrix(0,nrow=nrow(A),ncol=ncol(A))
    for(r in 1:(nrow(A) - 1)){ 
        # We will apply Householder to progressively the columns in A, decreasing 1 element at a time.
        x = A[r:nrow(A), r] 
        # We now get the vector v, starting with first entry = norm-2 of x[i] times 1
        # The sign is to avoid computational issues
        first = (sign(x[1]) * sqrt(sum(x^2))) +  x[1]
        # We get the rest of v, which is x unchanged, since e1 = [1, 0, 0, ..., 0]
        # We go the the last column / row, hence the if statement:
        v = if(length(x) > 1){c(first, x[2:length(x)])}else{v = c(first)}
        # Now we make the first entry unitary:
        w = v/first
        # Tau will be used in the Householder transform, so here it goes:
        t = as.numeric(t(w)%*%w) / 2
        # And the "reflectors" are stored as in the R qr()$qr function:
        reflectors[r: nrow(A), r] = w/t
        # The Householder tranformation is:
        I = diag(length(r:nrow(A)))
        H.transf = I - 1/t * (w %*% t(w))
        H_i  = diag(nrow(A))
        H_i[r:nrow(A),r:ncol(A)] = H.transf
        # And we apply the Householder reflection - we left multiply the entire A or Q
        A = H_i %*% A
        Q = H_i %*% Q
    }
    DECOMPOSITION = list("Q"= t(Q), "R"= round(A,7), 
            "compact Q as in qr()$qr"=  
            ((A*upper.tri(A,diag=T))+(reflectors*lower.tri(reflectors,diag=F))), 
            "reflectors" = reflectors,
            "rho"=c(apply(reflectors[,1:(ncol(reflectors)- 1)], 2, 
                function(x) sum(x^2) / 2), A[nrow(A),ncol(A)]))
    return(DECOMPOSITION)
}

Mari kita bandingkan ouput dengan fungsi bawaan R. Pertama fungsi buatan:

(H = House(X))
$Q
            [,1]        [,2]       [,3]       [,4]
[1,] -0.29329367 -0.73996967  0.5382474  0.2769719
[2,]  0.14829152 -0.65124800 -0.5656093 -0.4837063
[3,] -0.93923665  0.13835611 -0.1947321 -0.2465187
[4,] -0.09911084 -0.09580458 -0.5936794  0.7928072

$R
          [,1]       [,2]       [,3]      [,4]
[1,] -1.890006 -1.4517318  1.2524151 0.5562856
[2,]  0.000000 -0.9686105 -0.6449056 2.1735456
[3,]  0.000000  0.0000000 -0.8829916 0.5180361
[4,]  0.000000  0.0000000  0.0000000 0.4754876

$`compact Q as in qr()$qr`
            [,1]        [,2]       [,3]      [,4]
[1,] -1.89000649 -1.45173183  1.2524151 0.5562856
[2,] -0.14829152 -0.96861050 -0.6449056 2.1735456
[3,]  0.93923665 -0.67574886 -0.8829916 0.5180361
[4,]  0.09911084  0.03909742  0.6235799 0.4754876

$reflectors
            [,1]        [,2]      [,3] [,4]
[1,]  1.29329367  0.00000000 0.0000000    0
[2,] -0.14829152  1.73609434 0.0000000    0
[3,]  0.93923665 -0.67574886 1.7817597    0
[4,]  0.09911084  0.03909742 0.6235799    0

$rho
[1] 1.2932937 1.7360943 1.7817597 0.4754876

ke fungsi R:

qr.Q(qr(X))
            [,1]        [,2]       [,3]       [,4]
[1,] -0.29329367 -0.73996967  0.5382474  0.2769719
[2,]  0.14829152 -0.65124800 -0.5656093 -0.4837063
[3,] -0.93923665  0.13835611 -0.1947321 -0.2465187
[4,] -0.09911084 -0.09580458 -0.5936794  0.7928072

qr.R(qr(X))
          [,1]       [,2]       [,3]      [,4]
[1,] -1.890006 -1.4517318  1.2524151 0.5562856
[2,]  0.000000 -0.9686105 -0.6449056 2.1735456
[3,]  0.000000  0.0000000 -0.8829916 0.5180361
[4,]  0.000000  0.0000000  0.0000000 0.4754876

$qr
            [,1]        [,2]       [,3]      [,4]
[1,] -1.89000649 -1.45173183  1.2524151 0.5562856
[2,] -0.14829152 -0.96861050 -0.6449056 2.1735456
[3,]  0.93923665 -0.67574886 -0.8829916 0.5180361
[4,]  0.09911084  0.03909742  0.6235799 0.4754876

$qraux
[1] 1.2932937 1.7360943 1.7817597 0.4754876
Antoni Parellada
sumber
+1 tapi saya kira SPSS telah membangun fungsi untuk dekomposisi QR (atau apakah saya salah?), Jadi jika tujuannya adalah untuk membuat kode sesuatu di SPSS yang melibatkan QR, hampir tidak ada kebutuhan untuk mengimplementasikan algoritma QR secara manual.
Amuba kata Reinstate Monica
@amoeba Terima kasih. Karena kami sendirian, izinkan saya curhat: Penulis OP tidak memerlukan bantuan dari saya, tetapi saya mengambil komentarnya (di atas) tentang cara kerja bagian dalam (khususnya) fungsi R yang saya gunakan di pos pendamping sebagai tantangan pribadi yang menyenangkan, karena itu membuat saya menyadari betapa sedikitnya saya memahami implementasi faktorisasi QR ini yang dibangun ke dalam fungsi R. Karena sangat sulit bagi saya untuk menemukan referensi yang bagus, atau mendapatkan jawaban yang benar-benar membuat perbedaan dengan bertanya online, saya memposting apa yang saya dapatkan setelah lebih banyak upaya daripada yang ingin saya akui.
Antoni Parellada