Berechnung von Polynomkontrastvariablen

Bitte geben Sie mir eine Idee, wie eine kategoriale Variable (Faktor) effizient in die Menge der orthogonalen Polynomkontrastvariablen umcodiert werden kann.

Für viele Arten von Kontrastvariablen (z. B. Abweichung, einfach, Helmert usw.) lautet der Durchgang:

1. Stellen Sie die Kontrastkoeffizientenmatrix zusammen, die dem Typ entspricht.
2. Invers oder generalize-invers, um die Codematrix zu erhalten.

Beispielsweise:

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

Gibt es den gleichen oder einen ähnlichen Weg, um Polynomkontrastvariablen zu erhalten? Wenn ja, wie würde Matrix C aussehen und wie wird sie zusammengesetzt? Wenn nein , was nach wie vor ist die Art und Weise , das Polynom Kontrastvariablen effizient (zB durch Matrizenalgebra) zu berechnen.

Als Übergang zu meinem vorherigen Beitrag zu diesem Thema möchte ich einige vorläufige (wenn auch unvollständige) Untersuchungen der Funktionen hinter der linearen Algebra und verwandten R-Funktionen vorstellen. Dies soll in Arbeit sein.

Ein Teil der Undurchsichtigkeit der Funktionen hat mit der „Compact“ Form der Haus zu tun Zersetzung. Die Idee hinter der Householder-Zerlegung besteht darin, Vektoren über eine Hyperebene zu reflektieren, die durch einen Einheitsvektor u wie im folgenden Diagramm bestimmt wird, diese Ebene jedoch gezielt auszuwählen, um jeden Spaltenvektor der ursprünglichen Matrix A auf das e 1 zu projizieren Standardeinheitsvektor. Der normalisierte Norm-2 1 -Vektor u kann verwendet werden, um die verschiedenen Householder-Transformationen I - 2 zu berechnen$\mathbf {QR}$$\mathbf u$$\bf A$$\bf e_1$$1$$\bf u$ .$\mathbf{ I - 2\, uu^T x}$

Die resultierende Projektion kann ausgedrückt werden als

$\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}$

Der Vektor repräsentiert die Differenz zwischen den Spaltenvektoren x in der Matrix A , die wir zerlegen möchten, und den Vektoren y , die der durch u bestimmten Reflexion über den Unterraum oder "Spiegel" entsprechen .$\bf v$$\bf x$$\bf A$$\bf y$$\bf u$

$1$$\bf v$$\bf u$$\lVert u\rVert= 1$$1$$\bf w$ wir - weiterhin als Richtungsvektoren verwendet werden.

Das Schöne an der Methode ist, dass $\bf R$$\bf QR$$0$$\bf R$$\bf w$$1$$1$$\bf R$

qr()$qr$\mathbf R$ Matrix und der unteren dreieckigen "Speicher" -Matrix für die "modifizierten" Reflektoren verstanden werden.

$\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}$

Man würde annehmen, dass es in Ordnung wäre, diese Reflektoren unterhalb der Diagonale oder zu lagern$\bf w$$\bf R$$1$qr()$qr$\bf w$$\text{tau}$

können die Reflektoren alsReflektoren=ausgedrückt$\Large \tau = \frac{w^T\,w}{2}= \frac{\lVert \bf w \rVert}{2}$$\text{reflectors}=\large \bf w/\tau$$\bf R$$\bf QR$

$\bf w$$1$qr()qr()$qr$\tau$qr()$qraux$\rho=\frac{\sum \text{reflectors}^2}{2}= \frac{\bf w^Tw}{\tau^2} / 2$

$\bf w/\tau$

Der gesamte Code ist hier , aber da es bei dieser Antwort um den Schnittpunkt von Codierung und linearer Algebra geht, werde ich die Ausgabe zur Vereinfachung einfügen:

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

Jetzt habe ich die Funktion House()wie folgt geschrieben:

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

Vergleichen wir die Ausgabe mit den in R integrierten Funktionen. Zuerst die hausgemachte Funktion:

(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

zu den R-Funktionen:

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