Obliczanie wielomianowych zmiennych kontrastowych

11

Daj mi pomysł, jak efektywnie przekodować zmienną kategorialną (czynnik) do zestawu ortogonalnych wielomianowych zmiennych kontrastowych.

W przypadku wielu typów zmiennych kontrastu (np. Odchylenie, prosty, Helmert itp.) Przejście jest następujące:

  1. Skomponuj macierz współczynników kontrastu odpowiadającą typowi.
  2. Odwróć lub uogólnij-odwróć, aby uzyskać macierz kodów.

Na przykład:

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.)

Czy istnieje taki sam lub podobny sposób uzyskania wielomianowych zmiennych kontrastowych? Jeśli tak, to jak wyglądałaby matryca C i jak ją skomponować? Jeśli nie, co nadal jest sposobem na skuteczne obliczanie wielomianowych zmiennych kontrastowych (np. Za pomocą algebry macierzowej).

ttnphns
źródło
1
Patrzyłem na twoje pytanie po zweryfikowaniu (nawiasem mówiąc), że wyniki qr.qy()zgadzają się z ręcznymi obliczeniami, qr.Q(qr(X))a następnie Q%*%zw moim poście. Naprawdę zastanawiam się, czy mogę powiedzieć coś innego, aby odpowiedzieć na twoje pytanie bez powielania. Naprawdę nie chcę wykonywać złej roboty ... Przeczytałem wystarczająco dużo twoich postów, aby mieć dla ciebie duży szacunek ... Jeśli znajdę sposób na wyrażenie tej koncepcji bez kodu, tylko koncepcyjnie za pomocą algebry liniowej, Wrócę do tego. Cieszę się jednak, że odkryłem, że zgłębiłem zagadnienie o pewnej wartości. Najlepsze życzenia, Toni.
Antoni Parellada,
@Antoni, dziękuję. Moim celem jest możliwość kodowania (w SPSS, poprzez jego składnię). Czy można dowiedzieć się, jak działają wspomniane funkcje?
ttnphns

Odpowiedzi:

5

Jako fragment mojego wcześniejszego postu na ten temat chcę podzielić się pewną próbną (choć niepełną) eksploracją funkcji stojących za algebrą liniową i pokrewnymi funkcjami R. To ma być praca w toku.

Część nieprzejrzystości funkcji ma związek z „kompaktową” formą rozkładu Householder Ideą dekompozycji Householdera jest odbicie wektorów w poprzek hiperpłaszczyzny określonej przez wektor jednostkowy u, jak na poniższym schemacie, ale wybranie tej płaszczyzny w celowy sposób, aby rzutować każdy wektor kolumny oryginalnej macierzy A na e 1 standardowy wektor jednostki. Znormalizowany normą-2 1 wektor U mogą być wykorzystane do obliczenia RÓŻNYCH Transformations Householder I - 2QRuAe11u .I2uuTx

wprowadź opis zdjęcia tutaj

Wynikową projekcję można wyrazić jako

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

Wektor reprezentuje różnicę między wektorami kolumnowymi w macierzy , którą chcemy rozkładać, a wektorami odpowiadającymi odbiciu w podprzestrzeni lub „lustrze” wyznaczonemu przez .vxAUyu

Metoda zastosowana przez LAPACK uwalnia potrzebę przechowywania pierwszego wpisu w reflektorach Householder, zamieniając je w . Zamiast normalizować wektor do pomocą , jest to tylko pierwszy wpis konwertowany na ; jednak te nowe wektory - nazywaj je nadal mogą być używane jako wektory kierunkowe.v u u = 1 1 w1vuu=11w

Piękno tej metody jest to, że biorąc pod uwagę, że w sposób dekompozycja jest górna trójkątna, możemy rzeczywiście skorzystać z elementów poniżej przekątnej, aby wypełnić je z tymi reflektorów. Na szczęście wszystkie wiodące wpisy w tych wektorach są równe , co zapobiega problemowi z „sporną” przekątną macierzy: wiedząc, że wszystkie one są , nie trzeba ich uwzględniać i mogą dać przekątną wpisy do .Q R 0 R w 1 1 RRQR0Rw11R

Macierz „kompaktowego QR” w funkcji qr()$qrmożna z grubsza rozumieć jako dodatek macierzy i dolnej trójkątnej macierzy „przechowywania” dla „zmodyfikowanych” reflektorów.R

Projekcja Householder nadal będzie miała postać , ale nie będziemy pracować z ( ), ale raczej z wektorem , z czego tylko pierwszy wpis to , i u x = 1 w 1I2uuTxux=1w1

(1)I2uuTx=I2wwwTwx=I2wwTw2x .

Można by przypuszczać, że byłoby dobrze, aby przechowywać te reflektorów poniżej przekątnej lub wyłączeniem pierwszego wejścia , i nazywają to dzień. Jednak rzeczy nigdy nie są tak łatwe. Zamiast tego to, co jest przechowywane pod przekątną, to kombinacja i współczynników w transformacji Householdera wyrażonych jako (1), tak że definiując jako:wR1qr()$qrwtau

reflektory=w/τRQRτ=wTw2=w2 , reflektory można wyrazić jako . Te wektory „odblaskowe” są przechowywane bezpośrednio pod w tak zwanym „compact ”.reflectors=w/τRQR

Teraz jesteśmy o jeden stopień od wektorami, a pierwsza pozycja nie jest już , stąd wyjście będzie musiał zawierać klucz do ich przywrócenia ponieważ nalegają na wykluczeniu pierwszego wpisu z „odblaskowe” wektory do zmieści wszystko . Więc widzimy wartości na wyjściu? Cóż, nie, to byłoby przewidywalne. Zamiast tego w wynikach (gdzie ten klucz jest przechowywany) znajdujemy .1 τ ρ = reflektory 2w1qr()qr()$qrτqr()$qrauxρ=reflectors22=wTwτ2/2

Tak więc obramowane na czerwono poniżej, widzimy „reflektory” ( ), wyłączając ich pierwszy wpis.w/τ

wprowadź opis zdjęcia tutaj

Cały kod jest tutaj , ale ponieważ ta odpowiedź dotyczy przecięcia kodowania i algebry liniowej, dla ułatwienia wkleję wynik:


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

Teraz napisałem funkcję House()w następujący sposób:

   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)
}

Porównajmy wynik z wbudowanymi funkcjami R. Najpierw funkcja domowej roboty:

(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

do funkcji 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
źródło
+1, ale wydaje mi się, że SPSS ma wbudowane funkcje dekompozycji QR (czy się mylę?), Więc jeśli celem jest zakodowanie czegoś w SPSS, które obejmuje QR, nie ma potrzeby ręcznego wdrażania algorytmu QR.
ameba
@amoeba Dziękujemy. Ponieważ jesteśmy sami, pozwól, że ci się zwierzę: autor OP nie potrzebuje ode mnie żadnej pomocy, ale wziąłem jego komentarz (powyżej) na temat wewnętrznych funkcji (konkretnie) funkcji R, których użyłem w poście towarzyszącym jako osobiste wyzwanie zabawy, ponieważ uświadomiłem sobie, jak mało rozumiałem te implementacje faktoryzacji QR wbudowane w funkcje R. Ponieważ bardzo trudno było mi znaleźć jakieś dobre referencje lub uzyskać odpowiedzi, które naprawdę zmieniły się, pytając online, publikuję to, co otrzymałem po więcej wysiłku, niż chciałbym przyznać.
Antoni Parellada,