O que são polinômios ortogonais multivariados calculados em R?

12

Polinômios ortogonais em um conjunto univariado de pontos são polinômios que produzem valores nesses pontos de maneira que seu produto escalar e a correlação pareada sejam zero. R pode produzir polinômios ortogonais com a função poli .

A mesma função possui um polímero variante que produz polinômios ortogonais em um conjunto de pontos multivariados. De qualquer forma, os polinômios resultantes não são ortogonais no sentido de ter correlação zero em pares. De fato, como os polinômios de primeira ordem são apenas as variáveis ​​originais, os polinômios de primeira ordem não serão ortogonais, a menos que as variáveis ​​originais não estejam correlacionadas.

Então, minhas perguntas são:

  • Quais são os polinômios ortogonais multivariados calculados pelo polímero em R? Eles são apenas o produto dos polinômios ortogonais univariados? Para que são usados?
  • Podem existir polinômios ortogonais multivariados verdadeiros? Existe uma maneira fácil de produzi-los? Em R? Eles são realmente usados ​​em regressão?

Atualizar

Em resposta ao comentário de Superpronker, dou um exemplo do que quero dizer com polinômios não correlacionados:

> x<-rnorm(10000)
> cor(cbind(poly(x,degree=3)))
              1             2             3
1  1.000000e+00 -6.809725e-17  2.253577e-18
2 -6.809725e-17  1.000000e+00 -2.765115e-17
3  2.253577e-18 -2.765115e-17  1.000000e+00

A função poli retorna os polinômios ortogonais avaliados nos pontos x (aqui 10.000 pontos para cada polinômio). A correlação entre valores em diferentes polinômios é zero (com algum erro numérico).

Ao usar polinômios multivariados, as correlações são diferentes de zero:

> x<-rnorm(1000)
> y<-rnorm(1000)
> cor(cbind(polym(x,y,degree=2)))
              1.0           2.0           0.1         1.1           0.2
1.0  1.000000e+00  2.351107e-17  2.803716e-02 -0.02838553  3.802363e-02
2.0  2.351107e-17  1.000000e+00 -1.899282e-02  0.10336693 -8.205039e-04
0.1  2.803716e-02 -1.899282e-02  1.000000e+00  0.05426440  5.974827e-17
1.1 -2.838553e-02  1.033669e-01  5.426440e-02  1.00000000  8.415630e-02
0.2  3.802363e-02 -8.205039e-04  5.974827e-17  0.08415630  1.000000e+00

Portanto, não entendo em que sentido esses polinômios bivariados são ortogonais.

Atualização 2

Quero esclarecer o significado de "polinômios ortogonais" usados ​​na regressão, porque esse contexto pode ser de alguma maneira enganoso ao aplicar as idéias de polinômios ortogonais em intervalos conectados - como no comentário do Superpronker anterior.

Cito a Regressão Prática de Julian J. Faraway e Anova usando as páginas R 101 e 102:

Polinômios ortogonais contornam esse problema definindo

z1=a1+b1x
z2=a2+b2x+c2x2
z3=a3+b3x+c3x2+d3x3
etc. onde os coeficientes a, b, c ... são escolhidos para que ziT·zj=0quando ij . Os z são chamados polinômios ortogonais.

Por uma ligeira abuso de linguagem, aqui o autor utiliza zi tanto para o polinômio (como uma função) e para o vetor dos valores do polinômio leva nos pontos do conjunto x . Ou talvez nem seja um abuso de linguagem, porque desde o início do livro x tem sido o preditor (por exemplo, o conjunto de valores adotados pelo preditor).

Esse significado de polinômios ortogonais não é realmente diferente dos polinômios ortogonais em um intervalo. Podemos definir polinômios ortogonais da maneira usual (usando integrais) sobre qualquer conjunto mensurável com qualquer função de medida. Aqui temos um conjunto finito ( x ) e estamos usando produto de ponto em vez de integral, mas ainda são polinômios ortogonais se tomarmos nossa função de medida como o delta de Dirac nos pontos de nosso conjunto finito.

E em relação à correlação: produto pontual de vetores ortogonais em Rn (como a imagem de vetores ortogonais em um conjunto finito). Se o produto escalar de dois vetores for zero, a covariância é zero e se a covariância é zero, a correlação é zero. No contexto de modelos lineares é muito útil relacionar "ortogonal" e "não correlacionado", como no "desenho ortogonal de experimentos".

Pere
fonte
O que você quer dizer quando diz que os polinômios de um ponto não são correlacionados? Variáveis ​​estocásticas podem ser não correlacionadas; vetores podem ter um produto de ponto igual a zero.
Superpronker
Quando avaliados em um conjunto finito de pontos, obtemos um conjunto de valores para cada polinômio. Podemos calcular a correlação entre esses conjuntos de valores e, para os polinômios ortogonais, obtemos correlação zero. Como a correlação está relacionada à covariância e a covariância está relacionada ao produto de ponto, presumo que a correlação zero e o produto com ponto zero sejam equivalentes.
Pere
Desculpe se entendi mal, mas ainda não o entendo. A correlação é entre dois vetores em que você tem, digamos, N observações de cada um. Você quer dizer que o termo de primeira e segunda ordem não deve ser correlacionado? Então depende dos pontos em que você avalia. Em [-1; 1] eles não são, mas em [0; 1] eles são. Penso que a sua intuição para a relação entre ortogonalidade e falta de correlação não é precisa.
Superpronker
Atualizei a questão com isso, embora no contexto da regressão a ortogonalidade e a não correlação sejam quase sinônimos. Eu vinculei uma fonte. E sim, depende dos pontos que estamos avaliando. O primeiro argumento da ordem poli é o vetor de pontos que estamos avaliando e o primeiro passo dos meus exemplos é a geração de um vetor de pontos a avaliar. Na regressão, estamos interessados ​​em vetores ortogonais nos valores de nosso preditor.
Pere
Eu acho que o abuso da notação é mais problemático do que parece; a ortogonalidade de dois polinômios não é definida como o produto escalar sendo zero, não importa onde você avalie os polinômios. Pelo contrário, é que dois termos polinomiais (de ordens diferentes) devem ter produto com ponto zero no "sentido da função"; e produtos de ponto para funções são tipicamente integrais por alguma medida (por exemplo, função Peso). Veja en.m.wikipedia.org/wiki/Orthogonal_polynomials . Se eu estiver correto, isso explica a confusão. Mas no wiki há um comentário sobre a relação com os momentos.
Superpronker

Respostas:

5

Vamos explorar o que está acontecendo. Tenho certeza de que você já conhece a maior parte do material a seguir, mas para estabelecer notação e definições e tornar as idéias claras, abordarei o básico da regressão polinomial antes de responder à pergunta. Se desejar, Rvá para o cabeçalho "O que faz" cerca de dois terços do caminho nesta postagem e, em seguida, volte para todas as definições necessárias.

A configuração

Estamos considerando uma matriz modelo n×kX de potenciais variáveis ​​explicativas em algum tipo de regressão. Isso significa que estamos pensando nas colunas de X como sendo n vetores X1,X2,,Xk e estaremos formando combinações lineares delas, β1X1+β2X2++βkXk, para prever ou estimar uma resposta.

Às vezes, uma regressão pode ser aprimorada através da introdução de colunas adicionais criadas pela multiplicação de várias colunas de X , coeficiente por coeficiente. Esses produtos são chamados de "monômios" e podem ser escritos como

X1d1X2d2Xkdk

onde cada "potência" di é zero ou maior, representando quantas vezes cada X1 aparece no produto. Observe que X0 é um vetor n de coeficientes constantes ( 1 ) e X1=X si. Assim, monomios (como vectores) gerar um espaço de vector que inclui o espaço de coluna inicial de X. A possibilidade de que seja um espaço vetorial maior dá a esse procedimento um escopo maior para modelar a resposta com combinações lineares.

Pretendemos substituir a matriz original do modelo X por uma coleção de combinações lineares de monômios. Quando o grau de pelo menos um desses monômios excede 1, isso é chamado de regressão polinomial.

Gradings of polynomials

O grau de um monômio é a soma de seus poderes, d1+d2++dk. O grau de uma combinação linear de monômios (um "polinômio") é o maior grau entre os termos monomiais com coeficientes diferentes de zero. O grau tem um significado intrínseca, porque quando altera a base do espaço vectorial original, cada vector Xi é recém representado por uma combinação linear de todos os vectores; monômios X1d1X2d2Xkdktornam-se polinômios do mesmo grau; e consequentemente o grau de qualquer polinômio é inalterado.

O grau fornece uma "classificação" natural para essa álgebra polinomial: o espaço vetorial gerado por todas as combinações lineares de monômios em X de grau até e incluindo d+1, chamados de "polinômios de [ou até] grau d+1 em X, "estende-se do espaço vectorial de polinómios até grau d em X.

Usos da regressão polinomial

Freqüentemente, a regressão polinomial é exploratória no sentido de que não sabemos desde o início quais monômeros devem ser incluídos. O processo de criação de novas matrizes de modelo a partir de monômios e readequação da regressão pode precisar ser repetido várias vezes, talvez um número astronômico de vezes em algumas configurações de aprendizado de máquina.

Os principais problemas com essa abordagem são

  1. Os monômios geralmente introduzem quantidades problemáticas de "multicolinearidade" na nova matriz do modelo, principalmente porque os poderes de uma única variável tendem a ser altamente colineares. (A colinearidade entre os poderes de duas variáveis ​​diferentes é imprevisível, porque depende de como essas variáveis ​​estão relacionadas e, portanto, é menos previsível.)

  2. Alterar apenas uma coluna da matriz do modelo ou introduzir uma nova ou excluir uma pode exigir uma "reinicialização a frio" do procedimento de regressão, o que pode levar muito tempo para ser computado.

As classificações das álgebras polinomiais fornecem uma maneira de superar os dois problemas.

Polinômios ortogonais em uma variável

Dado um vetor de coluna única X, um conjunto de "polinômios ortogonais" para X é uma sequência de vetores de coluna p0(X),p1(X),p2(X), formados como combinações lineares de monômios somente em X - ou seja, como potências de X as seguintes propriedades:

  1. Para cada grau d=0,1,2,, os vetores p0(X),p1(X),,pd(X) geram o mesmo espaço vetorial que X0,X1,,Xd. (Observe que X0 é o vetor n de um e X1 é apenas X próprio).

  2. Os pi(X) são mutuamente ortogonais no sentido de que para ij,

    pi(X)pj(X)=0.

Normalmente, a matriz do modelo de substituição

P=(p0(X)p1(X)pd(X))
formada a partir desses monômios é escolhida para ser ortonormal , normalizando suas colunas para o comprimento da unidade:
PP=Id+1.
Porque o inverso de PP aparece na maioria das equações de regressão e o inverso da matriz de identidade Id+1 por si só, isso representa um enorme ganho computacional.

A ortonormalidade quase determina o pi(X). Você pode ver isso por construção:

  • O primeiro polinômio, p0(X), deve ser um múltiplo do n vetor 1=(1,1,,1) de comprimento da unidade. Existem apenas duas opções,±1/n1. É habitual escolher a raiz quadrada positiva.

  • O segundo polinômio, p1(X), deve ser ortogonal a 1. Ele pode ser obtido pela regressão X contra 1, cuja solução é o vector dos valores médios X = ˉ X 1 . Se os resíduos ε = X - X não são identicamente zero, eles dão apenas as duas soluções possíveis p 1 ( X ) = ± ( 1 / | | ε | |X^=X¯1.ϵ=XX^p1(X)=±(1/||ϵ||)ϵ.

...

  • Geralmente, pd+1(X) é obtido por regressão de Xd+1 contra p0(X),p1(X),,pd(X) e redimensionando os resíduos como um vetor de comprimento unitário. Existem duas opções de sinal quando os resíduos não são todos zero. Caso contrário, o processo termina: será infrutífera de olhar para quaisquer poderes mais altos de X. (Esse é um bom teorema, mas sua prova não precisa nos distrair aqui.)

Este é o processo de Gram-Schmidt aplicado à sequência intrínseca dos vetores X0,X1,,Xd,. Geralmente é calculado usando uma decomposição QR, que é quase a mesma coisa, mas calculada de maneira numericamente estável.

Essa construção gera uma sequência de colunas adicionais a serem consideradas, incluindo na matriz do modelo. A regressão polinomial em uma variável, portanto, geralmente prossegue adicionando elementos dessa sequência, um por um, em ordem, até que nenhuma melhoria adicional na regressão seja obtida. Como cada nova coluna é ortogonal às anteriores, inclusive ela não altera nenhuma das estimativas anteriores do coeficiente. Isso cria um procedimento eficiente e prontamente interpretável.

Polinômios em várias variáveis

A regressão exploratória (assim como o ajuste do modelo) geralmente prossegue considerando primeiro quais variáveis ​​(originais) incluir em um modelo; avaliando se essas variáveis ​​poderiam ser aumentadas, incluindo várias transformações delas, como monômios; e, em seguida, introduzindo "interações" formadas a partir de produtos dessas variáveis ​​e suas re-expressões.

A execução de tal programa começaria com a formação de polinômios ortogonais univariados nas colunas de X separadamente. Depois de selecionar um grau adequado para cada coluna, você introduziria interações.

Nesse ponto, partes do programa univariado são interrompidas. Que sequência de interações você aplicaria, uma a uma, até que um modelo adequado seja identificado? Além disso, agora que realmente entramos no campo da análise multivariável, o número de opções disponíveis e sua crescente complexidade sugerem que pode haver retornos decrescentes na construção de uma sequência de polinômios ortogonais multivariados. Se, no entanto, você tivesse essa sequência em mente, poderia calculá-la usando uma decomposição QR.


O que Rfaz

O software para regressão polinomial, portanto, tende a se concentrar na computação de sequências polinomiais ortogonais univariadas . É característico Restender esse suporte o mais automaticamente possível a grupos de polinômios univariados. É isso que polyfaz. (Seu companheiro polymé essencialmente o mesmo código, com menos sinos e assobios; as duas funções fazem as mesmas coisas.)

Especificamente, polycalculará uma sequência de polinômios ortogonais univariados quando receber um único vetor X, parando em um grau especificado d. (Se d é muito grande - e pode ser difícil prever quão grande é -, infelizmente, isso gera um erro.) Quando um conjunto de vetores X1,,Xk recebe um formato na matriz X, vai voltar

  1. Sequências de polinômios ortonormais p1(Xj),p2(Xj),,pd(Xj) para cada j com um grau máximo solicitado d. (Como o vetor constante p0(Xi) é comum a todas as variáveis ​​e é muito simples - geralmente é acomodado pela interceptação na regressão - Rnão se preocupa em incluí-lo.)

  2. Todas as interações entre os polinômios ortogonais até e incluindo os de grau d.

O passo (2) envolve várias sutilezas. Geralmente, por uma "interação" entre variáveis, queremos dizer "todos os produtos possíveis", mas alguns desses produtos possíveis terão graus maiores que d. Por exemplo, com 2 variáveis d=2, R calcula

p1(X1),p2(X1),p1(X2),p1(X1)p1(X2),p2(X2).

Rque não incluem as interacções maior grau p2(X1)p1(X2), p1(X1)p2(X2) (polinómios de grau 3) ou p1(X2)p2(X2) (um polinômio de grau 4). (Essa não é uma limitação séria, pois você pode facilmente calcular esses produtos ou especificá-los em um formulaobjeto de regressão .)

Outra sutileza é que nenhum tipo de normalização é aplicado a nenhum dos produtos multivariados. No exemplo, o único produto desse tipo é p1(X1)p1(X2). No entanto, não há garantia de que sua média seja zero e quase certamente não terá norma unitária. Neste sentido, é um verdadeiro "interacção" entre p1(X1) e p1(X2) e, como tal, pode ser interpretado como as interacções são geralmente em um modelo de regressão.

Um exemplo

Let's look at an example. I have randomly generated a matrix

X=(135624).
To make the calculations easier to follow, everything is rounded to two significant figures for display.

X1=(1,5,2)X10=(1,1,1)p0(X1)=(1,1,1)/3(0.58,0.58,0.58). The next step includes X11=X1 itself. To make it orthogonal to p0(X1), regress X1 against p0(X1) and set p1(X1) equal to the residuals of that regression, rescaled to have unit length. The result is the usual standardization of X1 obtained by recentering it and dividing by its standard deviation, p1(X1)=(0.57,0.79,0.23). Finally, X12=(1,25,4) is regressed against p0(X1) and p1(X1) and those residuals are rescaled to unit length. We cannot go any further because the powers of X1 cannot generate a vector space of more than n=3 dimensions. (We got this far because the minimal polynomial of the coefficients of X1, namely (t1)(t5)(t4), has degree 3, demonstrating that all monomials of degree 3 or larger are linear combinations of lower powers and those lower powers are linearly independent.)

The resulting matrix representing an orthonormal polynomial sequence for X1 is

P1=(0.580.570.590.580.790.200.580.230.78)

(to two significant figures).

In the same fashion, an orthonormal polynomial matrix for X2 is

P2=(0.580.620.530.580.770.270.580.150.80).

The interaction term is the product of the middle columns of these matrices, equal to (0.35,0.61,0.035). The full matrix created by poly or polym, then, is

P=(0.570.590.620.350.530.790.200.770.610.270.230.780.150.0350.80).

Notice the sequence in which the columns are laid out: the non-constant orthonormal polynomials for X1 are in columns 1 and 2 while those for X2 are in columns 3 and 5. Thus, the only orthogonality that is guaranteed in this output is between these two pairs of columns. This is reflected in the calculation of PP, which will have zeros in positions (1,2),(2,1),(3,5), and (5,3) (shown in red below), *but may be nonzero anywhere else, and will have ones in positions (1,1),(2,2),(3,3), and (5,5) (shown in blue below), but is likely not to have a one in the other diagonal positions ((4,4) in this example). Indeed,

PP=(1010.280.091010.0910.3110.09110.2500.280.30.250.50.320.091100.321).

When you inspect the P matrix shown in the question, and recognize that multiples of 1017 are really zeros, you will observe that this pattern of zeros in the red positions holds. This is the sense in which those bivariate polynomials are "orthogonal."

whuber
fonte
1
(+1) Great read as usual. I believe there is a small typo: You write that R calculates p1(X1)p2(X2) but shouldn't it be p1(X1)p1(X2)?
COOLSerdash
1
@Cool Good catch--fixed now.
whuber
1
Thanks for that great answer. The fact that the answer arrives a long after I lost hope in it being answered makes it a very enjoyable surprise.
Pere
And I think there is another small typo: I think "X1=X itself" in the 4th paragraph is intended to be "X1=X itself".
Pere
Completely right. I am grateful that you are reading the text so closely that you find these errors!
whuber