+3 Daumen
2,4k Aufrufe

Teil 2: https://www.mathelounge.de/529253/artikel-regression-berechnung-regr…

LinearRegression.png

1. Einführung

Mit der Methode der kleinsten Quadrate kann man zu einer gegebenen Menge von Punkten eine Funktion finden, die möglichst nahe an den Punkten verläuft. Wir werden in diesem Artikel eine allgemeine Formel zur Berechnung einer linearen Ausgleichsfunktion y=a0+a1xy=a_0+a_1x für endlich viele Messwerte herleiten, welche in ein Java-Programm gegossen und für diverse Probleme, bei denen Ausgleichsrechnungen eine Rolle spielen (physikalische Experimente, statistische Untersuchungen, ...), genutzt werden kann.

2. Datenpaare

Für unseren Algorithmus definieren wir einen neuen Datentyp „Data“, der Datenpaare repräsentiert. Bei Datenpaaren handelt es sich um 22-Tupel (geordnete Paare) der Form (x,y)(x,y) mit x,yRx,y\in\mathbb{R}. xx und yy werden als private, unveränderliche Objektvariablen vom Typ double gespeichert und im Konstruktor zugewiesen. Mit zwei Akzessormethoden (Getter) können diese Werte aus einem Data-Objekt gelesen werden. Eine mögliche Implementierung sieht wie folgt aus:

public class Data {
       
        // x-Wert
        private final double x;
       
        // y-Wert
        private final double y;
       
        // Konstruktor
        public Data(final double x, final double y){
                this.x = x;
                this.y = y;
        }
       
        // Getter für den x-Wert
        public double get_x(){
                return this.x;
        }
       
        // Getter für den y-Wert
        public double get_y(){
                return this.y;
        }
}

   Gegeben seien nn Datenpaare (x1,y1),(x2,y2),...,(xn,yn)(x_1,y_1),(x_2,y_2), ..., (x_n,y_n), die linear durch y=a0+a1xy=a_0+a_1x mit a0,a1Ra_0,a_1\in\mathbb{R} approximiert werden sollen. Gesucht ist die Lösung des linearen Gleichungssystems ATAx=ATbA^TAx=A^Tb. Dabei ist x=(a0a1)Tx= \left( \begin{matrix} a_0&a_1 \end{matrix} \right)^T , b=(y1y2...yn)Tb=\left( \begin{matrix} y_1&y_2&...&y_n \end{matrix} \right)^T und ARn×2A\in\mathbb{R}^{n\times 2} eine Matrix, die folgende Gestalt besitzt: A=(1x11x21xn) A=\left( \begin{matrix} 1&x_1\\ 1&x_2\\ \vdots&\vdots\\ 1&x_n \end{matrix} \right) Setzen wir A,AT,xA, A^T, x und bb in ATAx=ATbA^TAx=A^Tb, so erhalten wir: (11...1x1x2...xn)(1x11x21xn)ATA(a0a1)x=(11...1x1x2...xn)(y1y2yn)ATb \underbrace{\left( \begin{matrix} 1&1&...&1\\ x_1&x_2&...&x_n \end{matrix} \right)\cdot \left( \begin{matrix} 1&x_1\\ 1&x_2\\ \vdots&\vdots\\ 1&x_n \end{matrix} \right) }_{A^TA}\cdot \underbrace{\left( \begin{matrix} a_0\\ a_1 \end{matrix} \right)}_{x} = \underbrace{\left( \begin{matrix} 1&1&...&1\\ x_1&x_2&...&x_n\\ \end{matrix} \right)\cdot \left( \begin{matrix} y_1\\ y_2\\ \vdots\\ y_n \end{matrix} \right)}_{A^Tb} Wir berechnen ATAA^TA und erhalten: ATA=(11...1x1x2...xn)(1x11x21xn)=(nx1+x2+...+xnx1+x2+...+xnx12+x22+...+xn2)=(nk=1nxkk=1nxkk=1nxk2) A^TA=\left( \begin{matrix} 1&1&...&1\\ x_1&x_2&...&x_n \end{matrix} \right)\cdot \left( \begin{matrix} 1&x_1\\ 1&x_2\\ \vdots&\vdots\\ 1&x_n \end{matrix} \right) = \left( \begin{matrix} n & x_1+x_2+...+x_n\\ x_1+x_2+...+x_n&x^2_1+x^2_2+...+x^2_n \end{matrix} \right)\\ =\left( \begin{matrix} n & \sum\limits_{k=1}^{n}{x_k}\\ \sum\limits_{k=1}^{n}{x_k}&\sum\limits_{k=1}^{n}{x^2_k} \end{matrix} \right) Die Inverse (ATA)1\left(A^TA\right)^{-1}
der Matrix ATAA^TA, die zur Lösung des linearen Gleichungssystems genutzt wird, berechnen wir wie folgt. Sei φ : =det(ATA)1\varphi := \det(A^TA)^{-1} gegeben durch φ : =1det(ATA)=1nk=1nxk2(k=1nxk)2 \varphi := \dfrac{1}{\det\left(A^TA\right)} = \dfrac{1}{n\cdot \sum\limits_{k=1}^{n}{x^2_k} - \left(\sum\limits_{k=1}^{n}{x_k}\right)^2} Dann besitzt (ATA)1\left(A^TA\right)^{-1} folgende Gestalt: (ATA)1=φ(k=1nxk2k=1nxkk=1nxkn) \left(A^TA\right)^{-1}= \varphi\cdot \left( \begin{matrix} \sum\limits_{k=1}^{n}{x^2_k} & -\sum\limits_{k=1}^{n}{x_k}\\ -\sum\limits_{k=1}^{n}{x_k}&n \end{matrix} \right) Nun berechnen wir ATbA^Tb und erhalten: ATb=(11...1x1x2...xn)(y1y2yn)=(y1+y2+...+ynx1y1+x2y2+...+xnyn)=(k=1nykk=1nxkyk) A^Tb=\left( \begin{matrix} 1&1&...&1\\ x_1&x_2&...&x_n\\ \end{matrix} \right)\cdot \left( \begin{matrix} y_1\\ y_2\\ \vdots\\ y_n \end{matrix} \right)= \left( \begin{matrix} y_1+y_2+...+y_n\\ x_1\cdot y_1+x_2\cdot y_2+...+x_n\cdot y_n\\ \end{matrix} \right) = \left( \begin{matrix} \sum\limits_{k=1}^{n}{y_k}\\ \sum\limits_{k=1}^{n}{x_k\cdot y_k} \end{matrix} \right) Den Lösungsvektor x=(a0a1)x= \left(\begin{matrix} a_0\\ a_1 \end{matrix}\right) erhalten wir durch beidseitige Multiplikation der Inversen von ATAA^TA, denn (ATA)1ATA=Inx=(ATA)1ATbx=(ATA)1ATb\underbrace{\left(A^TA\right)^{-1}\cdot A^TA}_{=I_n} \cdot x = \left(A^TA\right)^{-1}A^Tb\Longleftrightarrow x=\left(A^TA\right)^{-1}A^Tb Der Lösungsvektor ist also: x=(a0a1)=φ(k=1nxk2k=1nxkk=1nxkn)(k=1nykk=1nxkyk)=φ((k=1nyk)(k=1nxk2)(k=1nxk)(k=1nxkyk)n(k=1nxkyk)(k=1nxk)(k=1nyk))x= \left(\begin{matrix} a_0\\ a_1 \end{matrix}\right)= \varphi\cdot \left( \begin{matrix} \sum\limits_{k=1}^{n}{x^2_k} & -\sum\limits_{k=1}^{n}{x_k}\\ -\sum\limits_{k=1}^{n}{x_k}&n \end{matrix} \right)\cdot \left( \begin{matrix} \sum\limits_{k=1}^{n}{y_k}\\ \sum\limits_{k=1}^{n}{x_k\cdot y_k} \end{matrix} \right) \\\\\\ =\varphi \cdot \left( \begin{matrix} \left(\sum\limits_{k=1}^{n}{y_k}\right)\cdot \left(\sum\limits_{k=1}^{n}{x^2_k}\right)-\left(\sum\limits_{k=1}^{n}{x_k}\right)\cdot \left(\sum\limits_{k=1}^{n}{x_k\cdot y_k}\right)\\ n\cdot \left(\sum\limits_{k=1}^{n}{x_k\cdot y_k}\right)-\left(\sum\limits_{k=1}^{n}{x_k}\right)\cdot \left(\sum\limits_{k=1}^{n}{y_k}\right) \end{matrix} \right)
Die Regressionskoeffizienten a0a_0 und a1a_1 der linearen Ausgleichsgerade y=a0+a1x y=a_0+a_1\cdot x können demnach direkt durch
a0=φ((k=1nyk)(k=1nxk2)(k=1nxk)(k=1nxkyk))a_0=\varphi\cdot \left(\left(\sum\limits_{k=1}^{n}{y_k}\right)\cdot \left(\sum\limits_{k=1}^{n}{x^2_k}\right)-\left(\sum\limits_{k=1}^{n}{x_k}\right)\cdot\left(\sum\limits_{k=1}^{n}{x_k\cdot y_k}\right)\right)
a1=φ(n(k=1nxkyk)(k=1nxk)(k=1nyk))a_1=\varphi\cdot\left(n\cdot \left(\sum\limits_{k=1}^{n}{x_k\cdot y_k}\right)-\left(\sum\limits_{k=1}^{n}{x_k}\right)\cdot \left(\sum\limits_{k=1}^{n}{y_k}\right)\right)
berechnet werden. Diese Funktion bietet die bestmöglichste Approximation an die gegebenen Messwerte mit dem kleinsten Betrag für den Fehlervektor ee nach der Methode der kleinsten Quadrate.

3. Implementierung

Wir stellen in einer beliebigen Klasse eine statische Utility-Funktion mit dem Namen „regression“, die als Argument eine Liste mit Data-Objekten erhält und die lineare Regression gemäß der des im vorangegangenen Abschnitt entwickelten Algorithmus vornimmt. Die Ergebnisparameter a0a_0 und a1a_1 werden in einem Array vom Typ double gespeichert und zurückgegeben.

public static double[] regression(final List<Data> data){
        final double[] result = new double[2];
       
        // Definition der Variablen für den Algorithmus
        double phi = 0;
        double sum_of_x_squared = 0;
        double sum_of_x = 0;
        double sum_of_y = 0;
        double sum_of_x_times_y = 0;
        double a_0 = 0;
        double a_1 = 0;
        final int m = data.size();
       
        // Berechnen der Summen
        for(final Data data_point : data){
                sum_of_x_squared += Math.pow(data_point.get_x(), 2);
                sum_of_x += data_point.get_x();
                sum_of_y += data_point.get_y();
                sum_of_x_times_y += (data_point.get_x() * data_point.get_y());
        }
               
        phi = 1.0/(m * sum_of_x_squared - Math.pow(sum_of_x, 2));
       
        // Berechnen der Regressionskoeffizienten
        a_0 = phi * (sum_of_y * sum_of_x_squared - sum_of_x * sum_of_x_times_y);
        a_1 = phi * (m * sum_of_x_times_y - sum_of_x * sum_of_y);
               
        // Ausgabe
        System.out.println(a_0 + (a_1 == 0 ? "" : (a_1 > 0 ? "+" : ""))+a_1+"*x");
       
        // Abspeichern der Ergebnisse
        result[0] = a_0;
        result[1] = a_1;
               
        return result;
}

4. Beispiel für die schrittweise Ermittlung der Regressionsgeraden

Gegeben seien die Messwerte {(x,y)x,yR} : ={(2,4),(1,2),(0,1),(1,0),(2,1)}\{(x,y)\mid x,y\in\mathbb{R}\} :=\left\{(-2,4),(-1,2),(0,1),(1,0),(2,1)\right\}. Wir suchen die lineare Ausgleichsgerade: y=a0x+a1y=a_0\cdot x+a_1 mit a0,a1Ra_0,a_1\in\mathbb{R} (die Form der Ausgleichsgeraden kann je nach Aufgabenstellung abweichen; beziehen Sie dies in Ihre Modellierung ein!). Mit den gegebenen Messwerten besitzt AA die Form A=(2111011121) bzw. AT=(2101211111)A=\left( \begin{matrix} -2&1\\ -1&1\\ 0&1\\ 1&1\\ 2&1\\ \end{matrix} \right)\text{ bzw. }A^T=\left( \begin{matrix} -2&-1&0&1&2\\ 1&1&1&1&1 \end{matrix} \right) Zu lösen ist folgendes lineare Gleichungssystem: AT(2111011121)ATA(a0a1)=AT(42101)ATb\underbrace{A^T\cdot \left( \begin{matrix} -2&1\\ -1&1\\ 0&1\\ 1&1\\ 2&1\\ \end{matrix} \right)}_{A^T A}\cdot \left(\begin{matrix} a_0\\ a_1 \end{matrix}\right) = \underbrace{A^T \cdot\left(\begin{matrix} 4\\ 2\\ 1\\ 0\\ 1 \end{matrix}\right) }_{A^T b} Wir berechnen ATAA^TA: ATA=(2101211111)T(2111011121)=(10005)A^TA=\left( \begin{matrix} -2&-1&0&1&2\\ 1&1&1&1&1 \end{matrix} \right)^T\cdot \left( \begin{matrix} -2&1\\ -1&1\\ 0&1\\ 1&1\\ 2&1\\ \end{matrix} \right)= \left(\begin{matrix} 10&0\\ 0&5 \end{matrix} \right) Anschließend ermitteln wir ATbA^T b: (2101211111)(42101)=(88) \left( \begin{matrix} -2&-1&0&1&2\\ 1&1&1&1&1 \end{matrix} \right) \cdot\left(\begin{matrix} 4\\ 2\\ 1\\ 0\\ 1 \end{matrix}\right)= \left(\begin{matrix} -8\\ 8 \end{matrix}\right) Wir lösen also folgendes lineare Gleichungssystem: (10005)(a0a1)=(88) \left(\begin{matrix} 10&0\\ 0&5 \end{matrix} \right) \cdot \left(\begin{matrix} a_0\\ a_1 \end{matrix}\right) = \left(\begin{matrix} -8\\ 8 \end{matrix}\right) Die Lösung kann mittels der Inversen oder dem Gauß-Jordan-Verfahren ermittelt werden. Man erhält schlussendlich: (a0a1)=(0.81.6)\left(\begin{matrix} a_0\\ a_1 \end{matrix}\right)= \left(\begin{matrix} -0.8\\ 1.6 \end{matrix}\right) Die Gleichung der Ausgleichsgerade lautet also y=0.8x+1.6=1.60.8xy=-0.8x+1.6=1.6-0.8x 
Im nächsten Teil werden wir dieses Ergebnis mit unserem Programm verifizieren.

geschlossen: Mathe-Artikel
Avatar von

Klasse Artikel!

unter  2. Datenpaare

Gegeben seien n Datenpaare

fällt mir die Matrixgleichung aweng zu plötzlich aus dem Text. Ich würde die Entwicklung, das Einsetzen der Punkte in die Ausgleichsgerade

 (xi,yi)g : y=a1  x+a0 \left(x_i, y_i \right) \in g: y \, = \, a_1 \; x + a_0

als (überbestimmtes) Gleichungssystem beschreiben

(a1  x1+a0=y1a2  x2+a0=y2...an  xn+a0=yn)\left(\begin{array}{r} a_1 \; x_1 + a_0 = y_1\\ a_2 \; x_2 + a_0 = y_2\\...\\ a_n \; x_n +a_0 = y_n\\\end{array}\right)

und explizit zeigen wie sich die Matrixgleichung A x = b begründet.

Warum jetzt die Matrixgleichung AT A x= AT b gerade eine Lösung des Problems (kleinste Quadrate) darstellt erschließt sich z.B.  wenn ich die Abweichungsquadrate (nomen est omen) berechne

Rk=(g(xk)yk)2R_k = \left(g \left(x_k \right) - y_k \right)^{2}

und in Summe

R(a1,a0) : =(a1  x1+a0y1)2+(a1  x2+a0y2)2+..+(a1  xn+a0yn)2R(a_1, a_0) \, := \, \left(a_1 \; x_1 + a_0 - y_1 \right)^{2} + \left(a_1 \; x_2 + a_0 - y_2 \right)^{2} +.. + \left(a_1 \; x_n + a_0 - y_n \right)^{2}

minimiere (partielle Ableitungen):

(dR(a1,a0)  daodR(a1,a0)  da1)=0\left(\begin{array}{r}dR \left(a_1, a_0 \right) \; da_o\\dR \left(a_1, a_0 \right) \; da_1\\\end{array}\right) = 0

was ein GLS ergibt

(n  a0+(x1+x2+..+xn)  a1=y1+y2+..+yn(1+x1+x2+..+xn)  a0+(x12+x22+xn2)  a1=x1  y1+x2  y2+..+xn  yn) \left(\begin{array}{r} n \; a_0 + \left( x_1 + x_2+..+x_n \right) \; a_1 &=& y_1 + y_2 +..+ y_n\\ \left(1 + x_1 + x_2 +..+ x_n \right) \;a_0 + \left(x_1^{2} + x_2^{2} + x_n^{2} \right) \; a_1 & =& x_1 \; y_1 + x_2 \; y_2 +..+ x_n \; y_n\\\end{array}\right)

das als Matrixgleichung geschrieben auf  AT A = AT b

(nx1+x2+..+xn1+x1+x2+..+xnx12+x22+..+xn2)(a0a1)=(y1+y2+..+ynx1y1+x2y2+..+xnyn)\left(\begin{array}{rr}n& x_1+x_2+..+x_n\\1+x_1+x_2+..+x_n& x_1^{2} + x_2^{2}+..+x_n^{2}\\\end{array}\right) \left(\begin{array}{r}a_0\\a_1\\\end{array}\right) = \left(\begin{array}{r} y_1+y_2+..+y_n\\ x_1 y_1 + x_2 y_2+..+ x_n y_n\\\end{array}\right)

führt und da treffen sich die Gedanken wieder....

https://ggbm.at/f4PUNWYZ

Tolle Ergänzung! Vielen Dank :-)

Ein anderes Problem?

Stell deine Frage