Berechnung von Polynomkontrastvariablen

11

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.

ttnphns
quelle
1
Ich habe mir Ihre Frage angesehen, nachdem ich (nebenbei) überprüft hatte, dass die Ergebnisse qr.qy()mit den manuellen Berechnungen übereinstimmen, die auf meinem Beitrag qr.Q(qr(X))gefolgt sind Q%*%z. Ich frage mich wirklich, ob ich etwas anderes sagen kann, um Ihre Frage ohne Doppelarbeit zu beantworten. Ich möchte wirklich keinen schlechten Job machen ... Ich habe genug von Ihren Posts gelesen, um viel Respekt vor Ihnen zu haben ... Wenn ich einen Weg finde, das Konzept ohne Code auszudrücken, nur konzeptionell durch lineare Algebra, Ich werde darauf zurückkommen. Ich bin jedoch froh, dass Sie meine Untersuchung des Themas von Wert gefunden haben. Beste Wünsche, Toni.
Antoni Parellada
@ Antoni, danke. Mein Ziel ist es, es codieren zu können (in SPSS anhand seiner Syntax). Ist es möglich herauszufinden, wie die von Ihnen erwähnten Funktionen funktionieren?
ttnphns

Antworten:

5

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 berechnenQRuAe11u .I2uuTx

Geben Sie hier die Bildbeschreibung ein

Die resultierende Projektion kann ausgedrückt werden als

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

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

1vuu=11w wir - weiterhin als Richtungsvektoren verwendet werden.

Das Schöne an der Methode ist, dass RQR0Rw11R

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

I2uuTxux=1w1

(1)I2uuTx=I2wwwTwx=I2wwTw2x

Man würde annehmen, dass es in Ordnung wäre, diese Reflektoren unterhalb der Diagonale oder zu lagernwR1qr()$qrwtau

können die Reflektoren alsReflektoren=ausgedrücktwerdenτ=wTw2=w2reflectors=w/τRQR

w1qr()qr()$qrτqr()$qrauxρ=reflectors22=wTwτ2/2

w/τ

Geben Sie hier die Bildbeschreibung ein

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
Antoni Parellada
quelle
+1, aber ich denke, dass SPSS Funktionen für die QR-Zerlegung eingebaut hat (oder irre ich mich?). Wenn also das Ziel darin besteht, etwas in SPSS zu codieren, das QR beinhaltet, muss der QR-Algorithmus kaum manuell implementiert werden.
Amöbe sagt Reinstate Monica
@amoeba Danke. Da wir alleine sind, möchte ich mich Ihnen anvertrauen: Der Autor des OP braucht keine Hilfe von mir, aber ich habe seinen Kommentar (oben) über das Innenleben (speziell) der R-Funktionen, die ich im Begleitbeitrag verwendet habe, aufgenommen als persönliche lustige Herausforderung, weil mir klar wurde, wie wenig ich diese in die R-Funktionen eingebauten Implementierungen der QR-Faktorisierung verstand. Da es für mich unglaublich schwierig war, gute Referenzen zu finden oder Antworten zu erhalten, die wirklich einen Unterschied machten, indem ich online fragte, poste ich, was ich nach mehr Mühe bekommen habe, als ich zugeben möchte.
Antoni Parellada