Cálculo de variáveis ​​de contraste polinomial

11

Por favor, me dê uma idéia de como recodificar com eficiência uma variável categórica (fator) no conjunto de variáveis ​​de contraste polinomial ortogonais.

Para muitos tipos de variáveis ​​de contraste (por exemplo, desvio, simples, Helmert etc.), o passe é:

  1. Componha a matriz dos coeficientes de contraste correspondente ao tipo.
  2. Inversa ou generalizada-inversa para obter a matriz de códigos.

Por exemplo:

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 a mesma maneira ou similar de obter variáveis ​​de contraste polinomial? Se sim, como seria a matriz C e como compor? Se não, qual ainda é o caminho para calcular as variáveis ​​de contraste polinomial de forma eficiente (por exemplo, álgebra matricial).

ttnphns
fonte
1
Eu olhei para sua pergunta depois de verificar (aliás) que os resultados de qr.qy()concordam com os cálculos manuais qr.Q(qr(X))seguidos por Q%*%zno meu post. Eu realmente me pergunto se posso dizer algo diferente para responder sua pergunta sem duplicação. Realmente não quero fazer um trabalho ruim ... Li suas postagens o suficiente para ter muito respeito por você ... Se eu encontrar uma maneira de expressar o conceito sem código, apenas conceitualmente através da álgebra linear, Eu voltarei a isso. Estou feliz, no entanto, que você tenha achado minha exploração da questão de algum valor. Muitas felicidades, Toni.
Antoni Parellada
@ Antoni, obrigado. Meu objetivo é poder codificá-lo (no SPSS, por sua sintaxe). É possível descobrir como funcionam as funções mencionadas?
precisa saber é o seguinte

Respostas:

5

Como segue o meu post anterior sobre esse tópico , quero compartilhar algumas tentativas (embora incompletas) de explorar as funções por trás da álgebra linear e funções R relacionadas. Isso deveria ser um trabalho em andamento.

Parte da opacidade das funções tem a ver com a forma "compacta" da decomposição Householder . A idéia por trás da decomposição do chefe de família é refletir vetores através de um hiperplano determinado por um vetor de unidade como no diagrama abaixo, mas escolher esse plano de maneira proposital para projetar todos os vetores de coluna da matriz original no vetor de unidade padrão . O vetor normalizado da norma-2 pode ser usado para calcular as diferentes transformações de Household .u A e 1 1 u I - 2QRuAe11uI2uuTx

insira a descrição da imagem aqui

A projeção resultante pode ser expressa como

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

O vetor representa a diferença entre os vetores de coluna na matriz que queremos decompor e os vetores correspondentes à reflexão no subespaço ou no "espelho" determinado por .vxAyu

O método usado pelo LAPACK libera a necessidade de armazenamento da primeira entrada nos refletores do agregado familiar, transformando-os em 's. Em vez de normalizar o vetor para com , é apenas a primeira entrada que é convertida em ; ainda assim, esses novos vetores - chamem-nos ainda podem ser usados ​​como vetores direcionais.1vuu=11w

A beleza do método é que, dado que em um decomposição é triangular superior, nós podemos realmente aproveitar os elementos em abaixo da diagonal para preenchê-los com estes refletores. Felizmente, as entradas principais nesses vetores são iguais a , impedindo um problema na diagonal "disputada" da matriz: sabendo que são todas as elas não precisam ser incluídas e podem render a diagonal para as entradas de .RQR0Rw11R

A matriz "QR compacta" na função qr()$qrpode ser entendida como aproximadamente a adição da matriz e a matriz triangular "de armazenamento" inferior para os refletores "modificados".R

A projeção do chefe de família ainda terá o formato , mas não trabalharemos com ( ), mas sim com um vetor , dos quais apenas a primeira entrada é garantida como eI2uuTxux=1w1

(1)I2uuTx=I2wwwTwx=I2wwTw2x .

Seria de supor que seria muito bem para armazenar estes refletores abaixo da diagonal ou excluindo a primeira entrada de , e chamá-lo um dia. No entanto, as coisas nunca são tão fáceis. Em vez disso, o que é armazenado abaixo da diagonal é uma combinação de e os coeficientes na transformação Household expressos como (1), de modo que, definindo como:wR1qr()$qrwtau

τ=wTw2=w2 , os refletores podem ser expressos como . Esses vetores "refletores" são aqueles armazenados logo abaixo de no chamado "compact ".reflectors=w/τRQR

Agora, estamos a um grau de distância dos vetores e a primeira entrada não é mais ; portanto, a saída de precisará incluir a chave para restaurá-los, pois insistimos em excluir a primeira entrada dos vetores "refletores" para encaixar tudo . Então, estamos vendo os valores na saída? Bem, não, isso seria previsível. Em vez disso, na saída de (onde essa chave está armazenada), encontramos .w1qr()qr()$qrτqr()$qrauxρ=reflectors22=wTwτ2/2

Assim, emoldurados em vermelho abaixo, vemos os "refletores" ( ), excluindo sua primeira entrada.w/τ

insira a descrição da imagem aqui

Todo o código está aqui , mas como essa resposta é sobre a interseção de codificação e álgebra linear, colarei a saída para facilitar:


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

Agora eu escrevi a função da House()seguinte maneira:

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

Vamos comparar a saída com as funções internas do R. Primeiro a função caseira:

(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

para as funções 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
fonte
+1, mas acho que o SPSS incorporou funções para decomposição do QR (ou estou errado?), Portanto, se o objetivo é codificar algo no SPSS que envolva o QR, dificilmente será necessário implementar manualmente o algoritmo QR.
Ameba diz Reinstate Monica
@amoeba Obrigado. Como estamos sozinhos, deixe-me confiar em você: O autor do OP não precisa de nenhuma ajuda minha, mas eu fiz o seu comentário (acima) sobre o funcionamento interno (especificamente) das funções R que usei no post complementar como um desafio divertido pessoal, porque me fez perceber o quão pouco entendi essas implementações da fatoração QR incorporadas às funções R. Como foi incrivelmente difícil para mim encontrar boas referências ou obter respostas que realmente fizeram a diferença perguntando on-line, estou postando o que obtive depois de mais esforço do que gostaria de admitir.
Antoni Parellada