Mehrfachregression

[F98]

Hallo Mathe-Fans,

folgendes Kennfeld:

X1    X2    X3    Y

1.2   3.1   2.    5.7
2.5   3.1   2.5   8.2
3.5   4.5   2.5   5.0
4.    4.5   3.    8.2
6.    5.    3.5   9.5

Berechnung des Lösungsvektors über ordinary least squares (OLS):

X  =

   [ 1.    1.2    3.1    2.  ] 
   [ 1.    2.5    3.1    2.5 ] 
   [ 1.    3.5    4.5    2.5 ] 
   [ 1.    4.     4.5    3.  ] 
   [ 1.    6.     5.     3.5 ]

y = [ 5.7, 8.2, 5.0, 8.2, 9.5 ]

b = inv(X' * X) * X' * y

b = [-2.1649851, -0.7144632, -1.7850398, 7.0941849]

daraus ergibt sich folgende lineare Näherungsfunktion:

y =-2.1649851 - 0.7144632x1 - 1.7850398x2 + 7.0941849x3

-> näherungsweise Fläche durch die Punktewolke (X1-X3) oben.

Solange also der Datenbestand rel. homogen ist, kann man gut genähert innerhalb der Flächenfunktion approximieren.
Das ganze sieht dann schon nicht mehr so schön aus, wenn die Fläche meinetwegen einen Sattel in der Mitte hat.

Ich suche nun eine Möglichkeit, wie man nichtlinear nähert, also eine Näherungsfunktion, die mir ein rel. genaue Abbild der "Landschaft" liefert.
U.U. polynomial, wobei man den Grad des Polynoms angeben kann und dabei eben mit mehreren Eingangsgrößen ein Kennfeld beschrieben kann.

Vorschläge?

MfG F98.

[otze]

Polynome eignen sich für sowas eigentlich nicht, da Polynome höheren Grades zu Überschwingern neigen.

Was für sowas öfter benutzt wird, sind Neuronale Netze die eine Art Verkettung von Sigmoid Funktionen darstellen. Das sind auch allgemeine funktionsapproximatoren wie Polynome, haben aber wesentlich bessere Eigenschaften.

Das Problem bei nichtlinearer optimierung ist allerdings, dass alles was komplexer als Polynome zweiten Grades ist, nen ziemlich heißes Eisen wird. Hauptproblem ist, dass es nicht mehr nur ein globales Optimum wie bei den linearen Verfahren gibt, sondern unter Umständen verdammt viele lokale Optima.

Was du also brauchst ist den richtigen Algorithmus, und da gibt es ziemlich viele, die alle gemeinsam haben, dass sie ekelhaft selbst zu implementieren sind, deswegen solltest du dir dann eine Bibliothek in der Programmiersprache deiner Wahl suchen.

Was sich anbieten würde, wäre zum beispiel RProp oder Verfahren die backpropagation of Error beherrschen(Quickprop). Als Alternative, wenn du mal eine richtig Ekelhafte Landschaft der Fehlerfunktion hast(viele lokale optima), kannst du noch evolutionäre Verfahren wie NEAT oder CMA-ES verwenden.
(Die Liste ist bei weitem nicht komplett. Aber das sind die Algorithmen die ich mal in der Hand hatte)

[F98]

Die Sache läuft also Deiner Meinung nach mehr auf evolutionäre Algorithmen raus?

Mir wäre es lieber, ich hätte ein Verfahren, was durch Matritzen-Operationen zum Ziel führt.

[otze]

Nein, ich bin nicht unbedingt für evolutionäre Algorithmen in dem Fall. Quickprop oder IProp sind hier das Mittel der Wahl - sie sind einfach um ein Vielfaches schneller. Evolutionäre Algorithmen sind die Holzhammer, die du dann verwendest, wenn die anderen Mittel versagen oder deine Fehlerfunktion nicht differenzierbar ist.

Wenn du nur reine Matrixoperationen auf den Daten durchführen willst, dann kommst du nicht weit. Matrizen bilden nur lineare Operationen ab, aber für die Problemklasse die du da beschreibst, brauchst du Operationen die eine starke Nichtlinearität einführen. Deswegen wirst du um sowas wie einen Gradientenabstieg auf der Fehlerfunktion nicht herumkommen(der natürlich lokal linear ist, weswegen dann doch oft wieder Matrizen zum Zug kommen, aber da das Ganze iterativ geschieht, kriegst du da dann doch wieder starke Nichtlinearität rein)

[F98]

Gibt es nicht auch einen Algorithmus, wo man z.B. sich über den Daumen gepeilt eine nichtlineare Gleichung vorgibt z.B.:

y = b1 * sin(x1) + b2 * sqrt(x1) + b3 * x2^2 + b4 * x2^3

und dann schrittweise die Koeffizienten genähert werden? Durch die Wichtung der einzelnen b könnte man dann sagen welches Teilglied rausfällt.

z.B.

b =

    0.2
    0.006
    -1.4
    1.8

Also die Zielgleichung dann so aussehen würde, da man meinetwegen b1 = 0.006 weglassen könnte, weil dessen Einfluss zu gering ist:

y = 0.2 * sin(x1) + 1.4 * x2^2 + 1.8 * x2^3

[knivil]

Least square error sollte das auch leisten koennen.

y =-2.1649851 * -0.7144632x1 * -1.7850398x2 * 7.0941849x3

Hast du hier * mit + verwechselt?

[F98]

Ja, habs korrigiert.

Hast Du ein konkretes Beispiel?

[knivil]

http://userpage.fu-berlin.de/~voelker/gnuplotkurs/gnuplotkurs.html unter "Polynom an Daten anpassen (least squares fit)" ist ein Beispiel fuer den eindimensionalen Fall. Wie es konkret im mehrdimensionalen ausschaut, kann ich jetzt aber nicht sagen.

[F98]

... leider nur für den eindimensionalen Fall und Polinomialregression. Ich benötige aber eine kompliziertere Variante.

[otze]

knivil schrieb:

Least square error sollte das auch leisten koennen.

jep, das ist die Grundlage der von mir genannten Algorithmen. Gnuplot macht dann am Ende auch nichts anderes als die Parameter mit einem der Standardalgorithmen zu berechnen.

F98 schrieb:

Gibt es nicht auch einen Algorithmus, wo man z.B. sich über den Daumen gepeilt eine nichtlineare Gleichung vorgibt [...] und dann schrittweise die Koeffizienten genähert werden?

Ja, und sie wurden bereits genannt. Die ganzen hier genannten Algorithmen tun nichts anderes. Gradientenabstieg ist nichts witeres als "schrittweise die Koeffizienten annähern".

Das Problem ist, dass die Gleichung ja den MSE optimieren muss:

E_{MSE}[f]=\sum_{i=1}^n (y\_i-f(x\_i))^2

Und der ist optimal, wenn die Ableitung von E[f]=0 ist. Leider ist das Finden der Nullstellen bei allem was komplexer als quadratische Funktionen ist eine ekelhafte Sache. Kannst ja gerne mal deine über den Daumen gepeilte nichtlineare Gleichung einsetzen, Ableiten und 0 setzen. Dann wirste schnell merken, dass das was du da machen willst eine echt komplizierte Angelegenheit ist(wenn du f(x)=ax^3+b setzt, kommt da am Ende ein Polynom 5. Grades raus für das du dann Nullstellen suchen darfst. Mit e-Funktionen wird das sofort noch ätzender. und der Sinus tötet dich komplett, weil er erstmal unendlich viele Nullstellen hinzufügt )

[knivil]

y = b1 * sin(x1) + b2 * sqrt(x1) + b3 * x2^2 + b4 * x2^3

Mir faellt grade auf, dass das auch nur ein lineares Problem ist, da die gesuchten Parameter b1, b2, b3 und b4 nur linear mit einfliessen. Mach es dir also nicht so schwer, linear ist einfach.

[otze]

knivil schrieb:

y = b1 * sin(x1) + b2 * sqrt(x1) + b3 * x2^2 + b4 * x2^3

Mir faellt grade auf, dass das auch nur ein lineares Problem ist, da die gesuchten Parameter b1, b2, b3 und b4 nur linear mit einfliessen. Mach es dir also nicht so schwer, linear ist einfach.

Wie wäre dann denn der Lösungsansatz wenn es linear wäre? Das Einzige was mir einfiele, wäre die Daten in die einzelnen Terme zu transformieren und für jeden Datensatz den neuen Datenpunkt (sin(x1),sqrt(x1),x2^2,x23)^T auszurechnen und dann in diesem 4Dimensionalen Raum linear zu optimieren. Bin mir aber grad nicht sicher, ob das wirklich das macht, was wir wollen.

Natürlich klappt das nicht, wenn wir zum Beispiel b1*sin(b2*x1) haben.

[F98]

Ok, wie wäre es denn, wenn man z.b. folgende Zerlegung vornimmt:

x1   x2   y
0.5  1.3  1.4
0.7  1.5  2.8
0.9  1.9  2.6  
1.2  2.3  3.7
1.5  4.2  4.3
Modellgleichung linearer Ansatz:

y = b0 + b1x1 + b2x2

und nun den linearen Ansatz "polinomial erweitert", man "mogelt" also beliebige neue Koeffizienten auf Basis von x1, x2 dazu:

y = b0 + b1x1 + b2x1² + b3x2 + b4x2²

Nun muss man für die neuen Koeffizienten eine Tabelle mit "neuen" Datenwerten auf Basis von x1, x2 ergänzen:

x1   x1²   x2   x2²    y
x1   x2    x3   x4
0.5  0.25  1.3  1.69   1.4
0.7  0.49  1.5  2.25   2.8
0.9  0.81  2.6  6.76   2.6
1.2  1.44  3.7  13.69  3.7
1.5  2.25  4.2  17.64  4.3

mit SciLab sieht das folgendermaßen aus:

x1 = [0.5,0.7,0.9,1.2,1.5]
x2 = [0.25,0.49,0.81,1.44,2.25]
x3 = [1.3,1.5,2.6,3.7,4.2]
x4 = [1.69,2.25,6.76,13.69,17.64]
y = [1.4,2.8,2.6,3.7,4.3]

X = [ones(5,1), x1', x2', x3', x4']

X  =

    1.    0.5    0.25    1.3    1.69   
    1.    0.7    0.49    1.5    2.25   
    1.    0.9    0.81    2.6    6.76   
    1.    1.2    1.44    3.7    13.69  
    1.    1.5    2.25    4.2    17.64

b = inv(X'*X)*X'*y'

b =

  - 1.1530301  
    16.200242  
  - 5.805805   
  - 3.9454127  
    0.6114774

X*b =

    1.4  
    2.8  
    2.6  
    3.7  
    4.3

y = -1.1530301 + 16.200242x1 - 5.805805x1² - 3.9454127x2 + 0.6114774x2²

Wäre doch ein schöner Ansatz einen ausreichenden Datenbestand (Kennfeld) mit einem Modellgleichung basierend auf Polynomen (o.ä.) mit beliebiger Genauigkeit, d.h. bel. neuen Koeffizienten, zu approximieren. Man könnte auch Wurzel und Sinusfunktionen usw. mit in die Gleichung einbeziehen. Letztendlich ist das alles dann nur eine Frage des Rechenaufwandes bezüglich der Matrizen.

[otze]

Wäre doch ein schöner Ansatz einen ausreichenden Datenbestand (Kennfeld) mit einem Modellgleichung basierend auf Polynomen (o.ä.) mit beliebiger Genauigkeit, d.h. bel. neuen Koeffizienten, zu approximieren. Man könnte auch Wurzel und Sinusfunktionen usw. mit in die Gleichung einbeziehen. Letztendlich ist das alles dann nur eine Frage des Rechenaufwandes bezüglich der Matrizen.

Ich hab das vorhin mal in meinen Vorlesungsunterlagen angeschaut. Das wird so auch zum Teil gemacht, aber so einfach, wie das klingt, ist es nicht.

Natürlich geht das alles, wenn man den Rechenaufwand außer betracht lässt - aber der ist Astronomisch.

um ein komplettes Polynom zweiten Grades mit 2D Eingabe abzubilden, brauchst du bereits folgende Basispolynome zur Linearisierung: (1,x,y,x^2,y2,x*y) für Polynome 3. Grades kannst du noch folgende Terme hinzufügen: (x^3,x2*y,x*y^2,y3). Linearisierung eines Polynoms dritten Grades benötigt also bereits einen 10 Dimensionalen Parameterraum. Und das alles bei 2D. Höhere Dimensionen der Eingabedaten werden richtig schlimm. Für ein Polynom von Grad k und eine n-dimensionale Eingabe lässt sich die Anzahl der Parameter so berechnen:

$P=\sum_{i=0}^k { n+k-1 \choose k }$

Und das ist exponentiell in n und k - weswegen die Technik bei vielen Anwendungen ausscheidet (man muss ja auch bedenken, dass danach noch Matrixmultiplikationen kommen, die dann richtig teuer werden...).

//edit nochn Problem: du brauchst dann natürlich auch mindestens P samples, sonst kannste die Matrix nicht invertieren.

[F98]

Dass das Ganze mit teurem Rechenaufwand verbunden ist, ist mir bewußt. Ich beschränke mich auf max. 4 Eingangsgrößen.

Ich nehme mal an, dass Du mit x, y das Selbe meinst, wie ich mit x1, x2, oder?

Warum müssen eigentlich Deiner Angabe nach die Terme x*y (x1*x2 ??) mit als Basispolynom in die Gleichung aufgenommen werden? Kann man das nicht auch einfach weglassen, so wie ich es in meinem Beispiel oben gemacht habe?

Könntest Du mir evtl. mal die Seiten aus Deinen Unterlagen per Mail zukommen lassen, um dass ich besser den math. Hintergrund verstehe?

[pasti/off]

http://stat.ethz.ch/teaching/lectures/FS_2009/CompStat

[otze]

F98 schrieb:

Warum müssen eigentlich Deiner Angabe nach die Terme x*y (x1*x2 ??) mit als Basispolynom in die Gleichung aufgenommen werden? Kann man das nicht auch einfach weglassen, so wie ich es in meinem Beispiel oben gemacht habe?

Nein, kannst du nicht. Ohne diese Parameter kannst du zum Beispiel bei quadratischen Polynomen die Funktion nicht mehr rotieren. Das ist schlecht, wenn deine Funktion eine Ellipsenförmige Parabel als Form hat, die schräg im Raum liegt - da kommst du einfach nicht hin ohne den x1*x2 Parameter. Stichwort: Kovarianzmatrix und Gaussverteilung. Da siehst du ziemlich gut, was diese Parameter leisten. Kannst du auch schnell selbst testen. Rotiere folgede Matrix um einen beliebigen Winkel:

A=[c1,0]
  [0,c2]

Schaue dir dann an, was mit den Parametern die vorher 0 waren passiert. das ist dein (x1*x2) Parameter.

Könntest Du mir evtl. mal die Seiten aus Deinen Unterlagen per Mail zukommen lassen, um dass ich besser den math. Hintergrund verstehe?

Uff. Das ist da nur kurz zusammengefasst gewesen weil das nur ein Randthema war um ein anderes Problem zu lösen. Ich glaube nicht, dass das so viel bringt.

[F98]

Ich hoffe, ich verstehe Dich nicht falsch.

Mein Ziel ist es nicht, die Matrix zu rotieren. Es geht mir nicht um die 3D-grafische Verarbeitung der Daten in einem Spiel, sondern es geht letzendlich nur um eine Wissensbasis, die bis dato in Form eines n-dimensionalen gerasterten, normierten Kennfeldes abgelegt ist. Eingangsgrößen sind x1-xn, Ausgangsgröße y. Zusammenhänge beschrieben als WENN X1, X2, X3 DANN Y.

Auf dieses Kennfeld werden nun die Eingangsdaten gegeben und zurück kommt die Ausgangsgröße (z.B. eine Stellgröße). Dahinter steht bis dato ein recht komplexer Algorithmus, der im n-dimensionalen Suchraum rödelt. Das kostet Zeit und v.a. Speicherplatz.

Einfacher wäre es oft, wenn man so ein Kennfeld, für technische Prozesse ausreichend, näherungsweise als Gleichung beschreibt. Das ist und muss kein Hexenwerk sein. Man setzt also die Eingangsgrößen einfach in die linearisierte Gleichung ein und bekommt (näherungsweise) y raus.

Schön ist es eben, wenn man nun nicht nur eine lineare Näherung als Fläche durch den Kennfeldraum hat, sondern auch einen Sattel oder leicht geschwungenen Ansstieg nachbilden kann, nichts kompliziertes wie zerklüftete 3D-Landschaften, denn techn. Prozesse sind rel. homogen.

Die Ermittlung der Koeffizienten der zugehörigen Gleichung kann da meinetwegen 5-10 min dauern ;), das macht man ja nur einmal im Vorfeld. Wenn man die Gleichung einmal hat, kann man diese rel. leicht auf einer C-Basierten Plattform verwenden. Sowas macht sich besser als einen komplexen Algo zu implementieren, der auch noch viel Speicher auf der Zielplattform verbraucht.

[otze]

F98 schrieb:

Mein Ziel ist es nicht, die Matrix zu rotieren. Es geht mir nicht um die 3D-grafische Verarbeitung der Daten in einem Spiel

Mir auch nicht. Offensichtlich habe ich es nicht geschafft, meine Intention rüber zu bringen.

plotte mal folgende Funktion:

x^T*A*x mit x=(x1,x2), A=die Matrix aus dem vorherigen post und c1!=c2

das lässt sich auflösen in:
x^T*A*x=(c1*x1, c2*x2)*x=c1*x1^2+c2*x22

und das ist, wie du zugeben musst eine ganz normale quadratische Funktion wie du sie auch hattest.

nun rotiere die Matrix A und löse die Gleichung nochmal auf. du wirst sehen, dass sie auf einmal komplexer wurde - der plot der Funktion aber bis auf die Rotation gleich ist. Das hat nichts mit 3D Verarbeitung zu tun, sondern dass man Funktionen im R^n entlang einer/mehrerer Achse(n) rotieren kann. Die Form bleibt die Selbe, nur eben rotiert. Warum ist das wichtig?
Stell dir vor, deine Daten sind primär so verteilt, dass A sie annähern kann - dann ist alles in Ordnung. Alternativ könnten sie aber auch genauso rotiert sein wie ichs hier gezeigt habe.
Dann kann aber A die Daten nicht mehr länger annähern, weil A keine Rotation kennt.

Nebenbei: So wie du dein Problem beschrieben hast, bin ich mehr denn je für neuronale Netze und entsprechende Optimierungsalgorithmen. Das ist die klassische nichtlineare Regressionsaufgabe und neuronale Netze sind da eben unangefochten State of the Art - warum mit weniger begnügen?

[F98]

otze schrieb:

Stell dir vor, deine Daten sind primär so verteilt, dass A sie annähern kann

Das tritt in der Regel ein, da die Daten rel. homogen verteilt sind und auch das Modell wenig/nicht zerklüftet (rel. ebene Fläche, Sattel)

otze schrieb:

Alternativ könnten sie aber auch genauso rotiert sein wie ichs hier gezeigt habe.

So wie ich das aber rauslesen kommt sowas gar nicht bei mir in Frage. Wenn eine Rotation der Daten im Raum um eine Achse stattfände müßte ja das Polynom komplett neu approximiert werden, da sich die Stutzstellen der Fläche gleichmäßig um eine Achse mit einem best. Drehwinkel verschoben hätten.

[otze]

So wie ich das aber rauslesen kommt sowas gar nicht bei mir in Frage. Wenn eine Rotation der Daten im Raum um eine Achse stattfände müßte ja das Polynom komplett neu approximiert werden, da sich die Stutzstellen der Fläche gleichmäßig um eine Achse mit einem best. Drehwinkel verschoben hätten.

Du hast immer noch nicht verstanden, was ich gemeint habe.

Mach erstmal das, was ich dir vorhin als beispiel gezeigt habe. splotte die Funktionen, schau sie dir an. Und dann überleg, ob der parameter sinnlos ist. Wenn er etwas macht, was deine Funktion ohne ihn nicht könnte, dann ist er es nicht.

so, ich suche mal direkt ein Bild raus:

http://www.weberconnect.de/bjt4.gif

ich hoffe, du weißt wie man das ließt. Kannst du diese Funktion ohne den x*y Parameter annähern? (das ist eine gedrehte quadratische Funktion)