Obliczanie wielomianowych zmiennych kontrastowych

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

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 - $\mathbf {QR}$$\mathbf u$$\bf A$$\bf e_1$$1$$\bf u$ .$\mathbf{ I - 2\, uu^T x}$

Wynikową projekcję można wyrazić jako

$\text{sign}(x_i=x_1)\times \lVert x \rVert \begin{bmatrix}1\\0\\0\\\vdots\\0\end{bmatrix}+\begin{bmatrix}x_1\\x_2\\x_3\\\vdots\\x_m\end{bmatrix}$

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 .$\bf v$$\bf x$$\bf A$$\bf y$$\bf u$

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 $1$$\bf v$$\bf u$$\lVert u\rVert= 1$$1$$\bf w$

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 $\bf R$$\bf QR$$0$$\bf R$$\bf w$$1$$1$$\bf R$

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.$\mathbf 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 $\mathbf{ I - 2\, uu^T x}$$\bf u$$\lVert \bf x \rVert=1$$\bf w$$1$

$\mathbf{ I - 2\, uu^T x}=\mathbf{ I - 2\, \frac{w}{\lVert w \rVert} \frac{w^T}{\lVert w \rVert} x}=\mathbf{ I - 2\, \frac{w\,w^T}{\lVert w \rVert^2}\,x}\tag{1}$ .

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:$\bf w$$\bf R$$1$qr()$qr$\bf w$$\text{tau}$

reflektory=w/τRQ$\Large \tau = \frac{w^T\,w}{2}= \frac{\lVert \bf w \rVert}{2}$ , reflektory można wyrazić jako . Te wektory „odblaskowe” są przechowywane bezpośrednio pod w tak zwanym „compact ”.$\text{reflectors}=\large \bf w/\tau$$\bf R$$\bf QR$

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 $\bf w$$1$qr()qr()$qr$\tau$qr()$qraux$\rho=\frac{\sum \text{reflectors}^2}{2}= \frac{\bf w^Tw}{\tau^2} / 2$

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

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