Was sind in R berechnete multivariate orthogonale Polynome?

12

Orthogonale Polynome in einer univariaten Menge von Punkten sind Polynome, die Werte auf diesen Punkten auf eine Weise erzeugen, dass ihr Punktprodukt und ihre paarweise Korrelation Null sind. R kann orthogonale Polynome mit der Funktion poly erzeugen .

Dieselbe Funktion hat ein Variantenpolym, das orthogonale Polynome auf einer multivariaten Punktmenge erzeugt. Auf jeden Fall sind die resultierenden Polynome nicht orthogonal im Sinne einer paarweisen Nullkorrelation. Tatsächlich sind Polynome erster Ordnung nur dann orthogonal, wenn die ursprünglichen Variablen nicht korreliert sind, da angenommen wird, dass sie die ursprünglichen Variablen sind.

Dann sind meine Fragen:

  • Was sind die durch polym in R berechneten multivariaten orthogonalen Polynome? Sind sie nur das Produkt der univariaten orthogonalen Polynome? Wofür werden sie benutzt?
  • Kann es echte multivariate orthogonale Polynome geben? Gibt es eine einfache Möglichkeit, sie herzustellen? In R? Werden sie tatsächlich in der Regression eingesetzt?

Aktualisieren

Als Antwort auf Superpronkers Kommentar gebe ich ein Beispiel dafür, was ich mit unkorrelierten Polynomen meine:

> 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

Die Poly-Funktion gibt die orthogonalen Polynome zurück, die in Punkten x ausgewertet werden (hier 10.000 Punkte für jedes Polynom). Die Korrelation zwischen Werten auf verschiedenen Polynomen ist Null (mit einigen numerischen Fehlern).

Bei Verwendung multivariater Polynome unterscheiden sich die Korrelationen von Null:

> 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

Daher verstehe ich nicht, in welchem ​​Sinne diese bivariaten Polynome orthogonal sind.

Update 2

Ich möchte die Bedeutung von "orthogonalen Polynomen", die in der Regression verwendet werden, klarstellen, da dieser Kontext irreführend sein kann, wenn die Ideen von orthogonalen Polynomen in zusammenhängenden Intervallen angewendet werden - wie im letzten Superpronker-Kommentar.

Ich zitiere Julian J. Faraways Practical Regression und Anova mit R Seiten 101 und 102:

Orthogonale Polynome umgehen dieses Problem, indem sie

z1=a1+b1x
z2=a2+b2x+c2x2
z3=a3+b3x+c3x2+d3x3
usw. wobei die Koeffizienten a, b, c ... so gewählt sind, dass ziT·zj=0wenn ij . Die z heißen orthogonale Polynome.

Bei leichtem Sprachmissbrauch verwendet der Autor hier zi sowohl für das Polynom (als Funktion) als auch für den Vektor der Werte, die das Polynom in den Punkten der Menge x annimmt . Oder vielleicht ist es gar kein Sprachmissbrauch, weil x seit Beginn des Buches der Prädiktor ist (z. B. die Menge der vom Prädiktor genommenen Werte).

Diese Bedeutung von orthogonalen Polynomen unterscheidet sich nicht von orthogonalen Polynomen in einem Intervall. Wir können orthogonale Polynome auf die übliche Weise (unter Verwendung von Integralen) über jede meßbare Menge mit jeder Meßfunktion definieren. Hier haben wir eine endliche Menge ( x ) und verwenden Punktprodukt anstelle von Integral, aber das sind immer noch orthogonale Polynome, wenn wir unsere Maßfunktion als das Dirac-Delta in den Punkten unserer endlichen Menge nehmen.

Und in Bezug auf die Korrelation: Punktprodukt von orthogonalen Vektoren in Rn (als das Bild eines orthogonalen Vektors auf einer endlichen Menge). Wenn das Punktprodukt zweier Vektoren Null ist, ist die Kovarianz Null, und wenn die Kovarianz Null ist, ist die Korrelation Null. Im Zusammenhang mit linearen Modellen ist es sehr nützlich, "orthogonal" und "unkorreliert" in Beziehung zu setzen, wie dies bei der "orthogonalen Versuchsplanung" der Fall ist.

Pere
quelle
Was meinen Sie, wenn Sie sagen, dass die Polynome zu einem bestimmten Zeitpunkt nicht korreliert sind? Stochastische Variablen können unkorreliert sein. Vektoren können ein Skalarprodukt von Null haben.
Superpronker
Bei einer Bewertung mit einer endlichen Menge von Punkten erhalten wir eine Menge von Werten für jedes Polynom. Wir können die Korrelation zwischen diesen Wertesätzen berechnen und für orthogonale Polynome erhalten wir eine Nullkorrelation. Da sich die Korrelation auf die Kovarianz und die Kovarianz auf das Skalarprodukt bezieht, gehe ich davon aus, dass die Korrelation Null und das Skalarprodukt Null gleichwertig sind.
Pere
Tut mir leid, wenn ich es falsch verstehe, aber ich folge immer noch nicht. Die Korrelation besteht zwischen zwei Vektoren, von denen Sie beispielsweise N Beobachtungen haben. Meinen Sie, dass der Term erster und zweiter Ordnung nicht korreliert sein sollte? Dann kommt es auf die Punkte an, an denen Sie bewerten. Auf [-1; 1] sind sie nicht, aber auf [0; 1] sind sie. Ich denke, Ihre Intuition für die Beziehung zwischen Orthogonalität und Unkorreliertheit ist nicht präzise.
Superpronker
Ich habe die Frage damit aktualisiert, obwohl Orthogonalität und Unkorreliertheit im Kontext der Regression beinahe Synonyme sind. Ich habe eine Quelle verlinkt. Und ja, es hängt von den Punkten ab, die wir bewerten. Das erste Argument der Ordnung poly ist der Vektor der Punkte, die wir bewerten, und der erste Schritt meiner Beispiele ist die Erzeugung eines zu bewertenden Punktevektors. Bei der Regression interessieren uns orthogonale Vektoren in den Werten unseres Prädiktors.
Pere
Ich denke, der Missbrauch der Notation ist problematischer, als es scheint; Die Orthogonalität von zwei Polynomen wird nicht als Punktprodukt definiert, das Null ist, unabhängig davon, wo Sie die Polynome auswerten. Es ist vielmehr so, dass zwei Polynomterme (unterschiedlicher Ordnung) ein Nullpunktprodukt im "Funktionssinn" haben sollten; und Punktprodukte für Funktionen sind typischerweise Integrale bezüglich eines Maßes (dh Gewichtsfunktion). Siehe en.m.wikipedia.org/wiki/Orthogonal_polynoms . Wenn ich richtig liege, erklärt dies die Verwirrung. Aber im Wiki gibt es einen Kommentar zum Verhältnis zu Momenten.
Superpronker

Antworten:

5

Lassen Sie uns untersuchen, was gerade passiert. Ich bin mir sicher, dass Sie den größten Teil des folgenden Materials bereits kennen, aber um Notationen und Definitionen festzulegen und die Ideen klarer zu machen, werde ich die Grundlagen der polynomialen Regression behandeln, bevor ich die Frage beantworte. Wenn Sie möchten, springen Sie zu der Überschrift "Was Rmacht", die ungefähr zwei Drittel des Weges in diesem Beitrag ausmacht, und überspringen Sie dann die Definitionen, die Sie möglicherweise benötigen.

Die Einstellung

Wir betrachten eine n×k Modellmatrix X potenzieller erklärender Variablen in einer Art Regression. Dies bedeutet, dass wir uns die Spalten von X als n Vektoren X1,X2,,Xk vorstellen und lineare Kombinationen daraus bilden, β1X1+β2X2++βkXk, um eine Antwort vorherzusagen oder zu schätzen.

Manchmal kann eine Regression verbessert werden, indem zusätzliche Spalten eingefügt werden, die durch Multiplizieren verschiedener Spalten von X mit Koeffizienten für Koeffizienten erstellt werden. Solche Produkte heißen "Monomials" und können gerne geschrieben werden

X1d1X2d2Xkdk

wobei jede "Leistung" di Null oder größer ist, was darstellt, wie oft jedes X1 im Produkt vorkommt. Man beachte, dass X0 ein n -Vektor konstanter Koeffizienten ( 1 ) ist und X1=X selbst. Somit erzeugen Monome (als Vektoren) einen Vektorraum, der den ursprünglichen Spaltenraum von X. Die Möglichkeit, dass es sich um einen größeren Vektorraum handelt, gibt diesem Verfahren einen größeren Spielraum für die Modellierung der Antwort mit linearen Kombinationen.

Wir beabsichtigen, die ursprüngliche Modellmatrix X durch eine Sammlung linearer Kombinationen von Monomen zu ersetzen . Wenn der Grad mindestens eines dieser Monome 1, übersteigt , spricht man von einer polynomialen Regression.

Abstufungen von Polynomen

Der Grad eines Monoms ist die Summe seiner Potenzen, d1+d2++dk. Der Grad einer linearen Kombination von Monomen (ein "Polynom") ist der größte Grad unter den Monomialtermen mit Koeffizienten ungleich Null. Der Grad hat eine intrinsische Bedeutung, denn wenn Sie die Basis des ursprünglichen Vektorraums ändern, wird jeder Vektor Xi durch eine lineare Kombination aller Vektoren neu dargestellt. Monome X1d1X2d2Xkdkdadurch werden Polynome des gleichen Grades; und folglich bleibt der Grad eines Polynoms unverändert.

Der Grad liefert eine natürliche "Abstufung" für diese Polynomalgebra: der Vektorraum, der durch alle linearen Kombinationen von Monomen in X mit einem Grad bis einschließlich d+1, die "Polynome mit einem Grad von [oder bis zu] d+1 in" genannt werden X, "erweitert den Vektorraum der Polynome bis zum Grad d in X.

Verwendung der polynomialen Regression

Polynom-Regression ist oft explorativ in dem Sinne, dass wir zu Beginn nicht wissen, welche Monome eingeschlossen werden sollen. Der Prozess des Erstellens neuer Modellmatrizen aus Monomen und des erneuten Anpassens der Regression muss möglicherweise mehrmals wiederholt werden, möglicherweise astronomisch oft in einigen Einstellungen für maschinelles Lernen.

Die Hauptprobleme bei diesem Ansatz sind

  1. Monome führen oft problematische Mengen an "Multikollinearität" in die neue Modellmatrix ein, hauptsächlich weil Potenzen einer einzelnen Variablen dazu neigen, hochkollinear zu sein. (Die Kollinearität zwischen Potenzen zweier verschiedener Variablen ist nicht vorhersehbar, da sie davon abhängt, wie diese Variablen zusammenhängen, und daher weniger vorhersehbar ist.)

  2. Das Ändern nur einer einzelnen Spalte der Modellmatrix oder das Einfügen oder Löschen einer neuen Spalte kann einen "Kaltstart" der Regressionsprozedur erfordern, dessen Berechnung möglicherweise viel Zeit in Anspruch nimmt.

Die Abstufungen von Polynomalgebren bieten eine Möglichkeit, beide Probleme zu überwinden.

Orthogonale Polynome in einer Variablen

Bei einem einzelnen Spaltenvektor X, eine Menge von "orthogonalen Polynomen" für X eine Folge von Spaltenvektoren p0(X),p1(X),p2(X), die allein in X als lineare Kombinationen von Monomen gebildet werden --das ist als Potenzen von X --mit den folgenden Eigenschaften:

  1. Für jeden Grad d=0,1,2,, die Vektoren p0(X),p1(X),,pd(X) erzeugen , um den gleichen Vektorraum als X0,X1,,Xd. (Beachten Sie, dass X0 der n Vektor von Einsen ist und X1 nur X selbst.)

  2. Das pi(X) sind zueinander orthogonal in dem Sinne , daß für ij,

    pi(X)pj(X)=0.

Üblicherweise wird die aus diesen Monomen gebildete Ersatzmodellmatrix

P=(p0(X)p1(X)pd(X))
orthonormal gewählt, indem ihre Spalten auf Einheitslänge normiert werden:
PP=Id+1.
Da die Inverse von PP erscheint in den meisten Regressionsgleichungen und die Inverse der Einheitsmatrix Id+1 Dies ist selbst ein riesiger Rechengewinn.

Die Orthonormalität bestimmt fast das pi(X). Sie können dies anhand der Konstruktion sehen:

  • Das erste Polynom p0(X), muss ein Vielfaches von dem betragen n -Vektor 1=(1,1,,1) der Einheitslänge. Es gibt nur zwei Möglichkeiten, ±1/n1. Es ist üblich, die positive Quadratwurzel zu wählen.

  • Das zweite Polynom, p1(X), muss orthogonal zu 1. Es kann durch Regression erhalten werden X gegen 1, deren Lösung der Vektor der Mittelwerte X = ˉ X 1 . Wenn die Residuen ε = X - X ist nicht gleich Null, sie geben den nur zwei möglichen Lösungen p 1 ( X ) = ± ( 1 / | | ε | |X^=X¯1.ϵ=XX^p1(X)=±(1/||ϵ||)ϵ.

...

  • Im Allgemeinen wird pd+1(X) erhalten, indem Xd+1 gegen p0(X),p1(X),,pd(X) und die Residuen neu skaliert werden, um einen Vektor der Längeneinheit zu ergeben. Es gibt zwei Vorzeichenoptionen, wenn die Residuen nicht alle Null sind. Andernfalls endet der Prozess: Es ist sinnlos, höhere Potenzen von X. (Dies ist ein netter Satz, aber sein Beweis braucht uns hier nicht abzulenken.)

Dies ist der Gram-Schmidt-Prozess , der auf die intrinsische Sequenz der Vektoren X0,X1,,Xd,. angewendet wird . Normalerweise wird es mit einer QR-Zerlegung berechnet , die sehr ähnlich ist, aber numerisch stabil berechnet wird.

Diese Konstruktion ergibt eine Folge zusätzlicher Spalten, die in die Modellmatrix aufgenommen werden müssen. Die polynomielle Regression in einer Variablen wird daher gewöhnlich fortgesetzt, indem Elemente dieser Sequenz nacheinander hinzugefügt werden, bis keine weitere Verbesserung der Regression mehr erzielt wird. Da jede neue Spalte zu den vorherigen orthogonal ist und auch keine der vorherigen Koeffizientenschätzungen ändert. Dies führt zu einer effizienten und leicht interpretierbaren Prozedur.

Polynome in mehreren Variablen

Die explorative Regression (sowie die Modellanpassung) erfolgt normalerweise, indem zunächst überlegt wird, welche (ursprünglichen) Variablen in ein Modell einbezogen werden sollen. Anschließend wird bewertet, ob diese Variablen durch Einbeziehen verschiedener Transformationen wie Monome erweitert werden könnten. und dann Einführen von "Wechselwirkungen", die aus Produkten dieser Variablen und ihren erneuten Ausdrücken gebildet werden.

Die Ausführung eines solchen Programms würde dann mit der getrennten Bildung univariater orthogonaler Polynome in den Spalten von X . Nachdem Sie für jede Spalte einen geeigneten Grad ausgewählt haben, führen Sie Interaktionen ein.

An diesem Punkt brechen Teile des univariaten Programms zusammen. Welche Abfolge von Interaktionen würden Sie nacheinander anwenden, bis ein geeignetes Modell gefunden ist? Darüber hinaus deuten die Anzahl der verfügbaren Optionen und ihre zunehmende Komplexität darauf hin, dass die Konstruktion einer Folge multivariater orthogonaler Polynome nachlässt, da wir uns nun wirklich mit der multivariablen Analyse befasst haben . Wenn Sie jedoch an eine solche Sequenz gedacht hätten, könnten Sie sie mithilfe einer QR-Zerlegung berechnen.


Was Rmacht

Software für die Polynomregression konzentriert sich daher in der Regel auf die Berechnung univariater orthogonaler Polynomsequenzen. Es ist charakteristisch R, diese Unterstützung so automatisch wie möglich auf Gruppen univariater Polynome auszudehnen. Das was polymacht. (Sein Begleiter polymist im Wesentlichen derselbe Code mit weniger Schnickschnack. Die beiden Funktionen führen die gleichen Aktionen aus.)

Insbesondere polywird eine Sequenz von univariaten orthogonalen Polynomen berechnet, wenn ein einzelner Vektor X, , der bei einem bestimmten Grad d. (Wenn d zu groß ist - und es kann schwierig sein, vorherzusagen, wie groß zu groß ist - wirft dies leider einen Fehler auf.) Wenn eine Menge von Vektoren X1,,Xk in Form einer Matrix X, es wird zurückkehren

  1. Sequenzen von orthonormalen Polynomen p1(Xj),p2(Xj),,pd(Xj) für jedes j bis zu einem angeforderten maximalen Grad d. (Da der konstante Vektor p0(Xi) allen Variablen gemeinsam ist und so einfach ist - er wird normalerweise vom Achsenabschnitt in der Regression Rberücksichtigt -, stört es nicht, ihn einzuschließen.)

  2. Alle Wechselwirkungen zwischen diesen orthogonalen Polynomen bis einschließlich derjenigen vom Grad d.

Schritt (2) beinhaltet mehrere Feinheiten. Normalerweise meinen wir mit "Wechselwirkung" zwischen Variablen "alle möglichen Produkte", aber einige dieser möglichen Produkte haben Grade größer als d. Beispielsweise mit 2 - Variablen und d=2, R berechnet

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

Rhabe nicht das höheren Grad Wechselwirkungen schließen p2(X1)p1(X2), p1(X1)p2(X2) (Polynome vom Grad 3) oder p1(X2)p2(X2) (ein Polynom 4. Grades). (Dies ist keine ernsthafte Einschränkung, da Sie diese Produkte problemlos selbst berechnen oder in einem Regressionsobjekt angeben können formula.)

Eine weitere feine Sache ist, dass auf keines der multivariaten Produkte eine Art Normalisierung angewendet wird. In dem Beispiel ist das einzige derartige Produkt p1(X1)p1(X2). Es gibt jedoch keine Garantie dafür, dass sein Mittelwert Null ist und es mit ziemlicher Sicherheit keine Einheitennorm gibt. In diesem Sinne handelt es sich um eine echte "Wechselwirkung" zwischen p1(X1) und p1(X2) als solche interpretiert werden kann, da sich Wechselwirkungen normalerweise in einem Regressionsmodell befinden.

Ein Beispiel

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.

The orthonormal polynomial sequence for the first column X1=(1,5,2) begins by normalizing X10=(1,1,1) to unit length, giving 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
quelle
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