Calcul des variables de contraste polynomiales

11

Veuillez me donner une idée de la façon de recoder efficacement une variable catégorielle (facteur) dans l'ensemble de variables de contraste polynomiales orthogonales.

Pour de nombreux types de variables de contraste (par exemple, déviation, simple, Helmert, etc.), la réussite est la suivante:

  1. Composez la matrice des coefficients de contraste correspondant au type.
  2. Inversez ou généralisez-l'inverse pour obtenir la matrice de codes.

Par exemple:

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

Existe-t-il la même méthode ou une méthode similaire pour obtenir des variables de contraste polynomiales? Si oui, à quoi ressemblerait la matrice C et comment la composer? Si non, quelle est encore la façon de calculer efficacement les variables de contraste polynomiales (par exemple par algèbre matricielle).

ttnphns
la source
1
J'ai regardé votre question après avoir vérifié (soit dit en passant) que les résultats qr.qy()sont en accord avec les calculs manuels qr.Q(qr(X))suivis de Q%*%zsur mon post. Je me demande vraiment si je peux dire autre chose pour répondre à votre question sans double emploi. Je ne veux vraiment pas faire un mauvais travail ... J'ai lu suffisamment de vos articles pour avoir beaucoup de respect pour vous ... Si je trouve un moyen d'exprimer le concept sans code, juste conceptuellement à travers l'algèbre linéaire, J'y reviendrai. Je suis heureux, cependant, que vous ayez trouvé mon exploration de la question d'une certaine valeur. Meilleurs voeux, Toni.
Antoni Parellada
@Antoni, merci. Mon objectif est de pouvoir le coder (en SPSS, par sa syntaxe). Est-il possible de comprendre comment fonctionnent les fonctions que vous mentionnez?
ttnphns

Réponses:

5

Comme enchaînement à mon post précédent sur ce sujet, je veux partager une exploration provisoire (quoique incomplète) des fonctions derrière l'algèbre linéaire et les fonctions R connexes. C'est censé être un travail en cours.

Une partie de l'opacité des fonctions est liée à la forme "compacte" de la décomposition Householder . L'idée derrière la décomposition du Householder est de refléter les vecteurs à travers un hyperplan déterminé par un vecteur-unité comme dans le diagramme ci-dessous, mais en choisissant ce plan de manière ciblée afin de projeter chaque vecteur de colonne de la matrice d'origine sur le vecteur d'unité standard . Le vecteur normalisé-2 normalisé peut être utilisé pour calculer les différentes transformations de Householder .u A e 1 1 u I - 2QRuAe11uI2uuTx

entrez la description de l'image ici

La projection résultante peut être exprimée comme

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

Le vecteur représente la différence entre les vecteurs colonne dans la matrice que nous voulons décomposer et les vecteurs correspondant à la réflexion à travers le sous-espace ou "miroir" déterminée par .vxAyu

La méthode utilisée par LAPACK libère le besoin de stockage de la première entrée dans les réflecteurs Householder en les transformant en . Au lieu de normaliser le vecteur en avec , c'est juste la première entrée qui est convertie en ; pourtant, ces nouveaux vecteurs - appelez-les peuvent toujours être utilisés comme vecteurs directionnels.1vuu=11w

La beauté de la méthode est que, étant donné que dans une décomposition est triangulaire supérieur, nous pouvons réellement profiter des éléments dans sous la diagonale pour les remplir avec ces réflecteurs . Heureusement, les entrées principales de ces vecteurs sont toutes égales à , ce qui évite un problème dans la diagonale "contestée" de la matrice: sachant qu'elles sont toutes elles n'ont pas besoin d'être incluses et peuvent donner la diagonale aux entrées de .RQR0Rw11R

La matrice "QR compacte" dans la fonction qr()$qrpeut être comprise en gros comme l'addition de la matrice et de la matrice triangulaire "de stockage" inférieure pour les réflecteurs "modifiés".R

La projection Householder aura toujours la forme , mais nous ne travaillerons pas avec ( ), mais plutôt avec un vecteur , dont seule la première entrée est garantie à , etI2uuTxux=1w1

(1)I2uuTx=I2wwwTwx=I2wwTw2x .

On pourrait supposer qu'il serait très bien de stocker ces réflecteurs sous la diagonale ou excluant la première entrée de , et l'appeler un jour. Cependant, les choses ne sont jamais aussi faciles. Au lieu de cela, ce qui est stocké en dessous de la diagonale dans est une combinaison de et des coefficients dans la transformation Householder exprimée comme (1), de sorte que, définissant comme:wR1qr()$qrwtau

τ=wTw2=w2 , les réflecteurs peuvent être exprimés sous la forme . Ces vecteurs "réflecteurs" sont ceux stockés juste sous dans le soi-disant "compact ".reflectors=w/τRQR

Maintenant, nous sommes à un degré des vecteurs , et la première entrée n'est plus , d'où la sortie de devra inclure la clé pour les restaurer car nous insistons pour exclure la première entrée des vecteurs "réflecteur" à s'adapter à tout . Voyons-nous donc les valeurs dans la sortie? Eh bien, non, ce serait prévisible. Au lieu de cela, dans la sortie de (où cette clé est stockée), nous trouvons .w1qr()qr()$qrτqr()$qrauxρ=reflectors22=wTwτ2/2

Ainsi encadrés en rouge ci-dessous, on voit les "réflecteurs" ( ), à l'exclusion de leur première entrée.w/τ

entrez la description de l'image ici

Tout le code est ici , mais comme cette réponse concerne l'intersection du codage et de l'algèbre linéaire, je vais coller la sortie pour plus de facilité:


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

Maintenant, j'ai écrit la fonction House()comme suit:

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

Comparons la sortie aux fonctions intégrées R. D'abord la fonction maison:

(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

aux fonctions 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
la source
+1 mais je suppose que SPSS a intégré des fonctions de décomposition QR (ou je me trompe?), Donc si le but est de coder quelque chose dans SPSS qui implique QR, il n'est guère nécessaire d'implémenter manuellement l'algorithme QR.
amibe dit Réintégrer Monica
@amoeba Merci. Puisque nous sommes seuls, permettez-moi de vous confier: l'auteur de l'OP n'a besoin d'aucune aide de ma part, mais j'ai pris son commentaire (ci-dessus) sur le fonctionnement interne (spécifiquement) des fonctions R que j'ai utilisées dans le message d'accompagnement comme un défi amusant personnel, car cela m'a fait réaliser à quel point je comprenais mal ces implémentations de la factorisation QR intégrées aux fonctions R. Comme il m'a été extrêmement difficile de trouver de bonnes références ou d'obtenir des réponses qui ont vraiment fait la différence en demandant en ligne, je publie ce que j'ai obtenu après plus d'efforts que je ne veux l'admettre.
Antoni Parellada