Wer ist online?
Buchempfehlung
FreeBASIC-Chat
FreeBASIC-Nachrichten jetzt auch über Twitter erhalten. Follow us!Code-Beispiel
X/Y-Regressionsanalyse
| Lizenz: | Erster Autor: | Letzte Bearbeitung: |
| k. A. |  TJF |  08.11.2011 |
Dieses Beispiel betrifft das mathematische Verfahren der 
Regressionsanalyse von Messdaten (auch Anpassung, Parameterschätzung, Ausgleichsrechnung oder Fitting genannt).
Aus einer Anzahl von diskreten Messwerten y_i an den Stellen x_i wird eine Fitting-Funktion y = F(x) berechnet, welche für weitere mathematische Bearbeitungen verwendet werden kann, z.B. zur Interpolation in Bereichen ohne Messdaten oder auch zur Bildung von Ableitungen der Messdaten. Übliche Verfahren betreffen entweder die Abhängigkeit von nur einer Eingangsgröße (y = F(x)) oder wenige Verfahren betreffen mehrere Eingangsgrößen (y = F(x1, x2, ..., xn) mit ausschließlich linearer Abhängigkeit der Eingangsgrößen, vgl. auch fbmath - Kap. 14.
Das hier vorgestellte Verfahren verwendet zwei Eingangsgrößen (x1, x2) und erzeugt eine Fitting-Funktion y = F(x1, x2) mit beliebiger Ordnung (= Potenz der Ansatzfunktionen, beliebig im Rahmen der Rechengenauigkeit).
Als Anwendungsbeispiel sei auf die Karten der Wettervorhersage verwiesen. Dort sind Linien eingezeichnet, die Hoch- und Tiefdruckgebiete kennzeichnen. Diese Isobaren werden aus Messwerten ermittelt, welche an diskreten Orten bestimmt werden. Es sind also Druckmesswerte P an den Orten (X, Y) verfügbar, aus denen dann der Verlauf der Isobaren ermittelt wird, also P(X, Y) = F(x1, x2) = y.
Alle Ausgleichsrechnungen lassen sich durch Variationsprobleme beschreiben, bei denen die Gewichtung von Ansatzfunktionen durch Parameter so gewählt wird, dass die Abweichung von den Messwerten minimal wird. Im Falle von mehreren Eingangsgrößen ist es jedoch nicht möglich, beliebige Ansatzfunktionen zu wählen, da sonst ggf. das Variationsproblem singulär wird, also unlösbar ist.
In diesem Code-Beispiel wird ein Produktansatz von Legendre Polynomen verwendet, welche im Intervall -1 bis 1 orthogonal sind. Daher ist bei hinreichender Verteilung der Messdaten das Variationsproblem immer lösbar. (Jedoch ist eine Koordinatentransformation unvermeidbar und die Fitting-Funktion sollte nur in dem zuvor definierten Gültigkeitsbereich angewendet werden.)
Zunächst der Quelltext für die Struktur XYRegression. Er verwendet das Code-Beispiel Legendre Polynome, welches wiederum libFBla verwendet. Der Quelltext sollte unter dem Namen "XYRegression.bas" gespeichert werden.
' This is file "XYRegression.bas"
' (C) LGPLv2 by Thomas[ dot ]Freiherr{ at }gmx{ dot }net
#INCLUDE ONCE "LegendrePolynom.bas"
#DEFINE XYR_INDEX(_L_,_H_) _L_ += 1 : IF _L_ >= Basis THEN _L_ = 0 : _H_ += 1
TYPE XYRegression
AS STRING Err_
AS LA_S StdDev ' standard deviation
DECLARE CONSTRUCTOR(BYVAL O AS UINTEGER, BYVAL V AS LA_M PTR, _
BYVAL ChX1 AS UINTEGER = 0, _
BYVAL ChX2 AS UINTEGER = 1, _
BYVAL ChY AS UINTEGER = 2, _
BYVAL StAbw AS UINTEGER = 0, _
BYVAL Xn AS LA_S = 0.0, BYVAL Xx AS LA_S = 0.0, _
BYVAL Yn AS LA_S = 0.0, BYVAL Yx AS LA_S = 0.0)
DECLARE DESTRUCTOR()
DECLARE FUNCTION Val_(BYVAL X AS LA_S, BYVAL Y AS LA_S) AS LA_S
DECLARE FUNCTION Abw_(BYVAL V AS LA_M PTR, _
BYVAL ChX1 AS UINTEGER = 0, _
BYVAL ChX2 AS UINTEGER = 1, _
BYVAL ChY AS UINTEGER = 2) AS LA_S
Private:
AS LegendrePolynom PTR Leg ' polynoms to use
AS LA_V Res ' vector to store parameters
AS LA_S Xmin, Ymin, Dx, Dy ' for transformation interval [-1, 1]
AS UINTEGER Basis, Lines, C1, C2, Cy ' for indexing
DECLARE SUB Fit_(BYVAL V AS LA_M PTR, BYVAL Mo AS INTEGER = 0)
END TYPE
' computes the Regression of Order from values V
' ChX1, ChX2 = columns of X values (defaults to 0, 1)
' ChY = column of Y values (defaults to 2)
' StAbw: if <> 0 then calculate the standard deviation
' Xn, Xm: minimum and maximum X1 values
' Yn, Ym: minimum and maximum X2 values (if Yn = Ym -> find extrema in V)
CONSTRUCTOR XYRegression(BYVAL Order AS UINTEGER, BYVAL V AS LA_M PTR, _
BYVAL ChX1 AS UINTEGER = 0, _
BYVAL ChX2 AS UINTEGER = 1, _
BYVAL ChY AS UINTEGER = 2, _
BYVAL StAbw AS UINTEGER = 0, _
BYVAL Xn AS LA_S = 0.0, BYVAL Xx AS LA_S = 0.0, _
BYVAL Yn AS LA_S = 0.0, BYVAL Yx AS LA_S = 0.0)
Basis = Order + 1
C1 = ChX1
C2 = ChX2
Cy = ChY
VAR az = V->Row_() + 1
Lines = Basis ^ 2
IF az < Lines THEN Err_ = "not enough values" : EXIT CONSTRUCTOR
VAR i = V->Col_()
IF i < C1 ORELSE _
i < C2 ORELSE _
i < Cy THEN Err_ = "column missmatch" : EXIT CONSTRUCTOR
IF Yn = Yx THEN ' search extrema
Xmin = V->Val_(0, C1) : Ymin = V->Val_(0, C2)
VAR Xmax = Xmin, Ymax = Ymin
FOR i AS INTEGER = 1 TO V->Row_
IF V->Val_(i, C1) < Xmin THEN Xmin = V->Val_(i, C1)
IF V->Val_(i, C1) > Xmax THEN Xmax = V->Val_(i, C1)
IF V->Val_(i, C2) < Ymin THEN Ymin = V->Val_(i, C2)
IF V->Val_(i, C2) > Ymax THEN Ymax = V->Val_(i, C2)
NEXT
Dx = Xmax - Xmin : IF Dx THEN Dx = 2.0 / Dx
Dy = Ymax - Ymin : IF Dy THEN Dy = 2.0 / Dy
ELSE ' compute transformation
Xmin = Xn : Ymin = Yn
Dx = Xx - Xn : IF Dx THEN Dx = 2 / Dx
Dy = Yx - Yn : IF Dy THEN Dy = 2 / Dy
END IF
IF Dx = 0 THEN Err_ = "no X difference" : EXIT CONSTRUCTOR
IF Dy = 0 THEN Err_ = "no Y difference" : EXIT CONSTRUCTOR
Leg = NEW LegendrePolynom(Basis - 1)
VAR w = LA_M(Lines, Lines), r = LA_V(Lines)
Lines -= 1
VAR x = LA_V(az), y = LA_V(az)
az -= 1 : Xmin += 1 / Dx : Ymin += 1 / Dy
FOR i = 0 TO az ' transform into [-1, 1]
x.Val_(i) = (V->Val_(i, C1) - Xmin) * Dx
y.Val_(i) = (V->Val_(i, C2) - Ymin) * Dy
NEXT
VAR m = 0, n = 0
DIM AS LA_S sum
FOR i = 0 TO Lines ' build linear equation system
VAR j = 0, k = 0, p = 0, q = 0
FOR k = 0 TO Lines ' left side
sum = 0.0
FOR j = 0 TO az
sum += Leg->Val_(x.Val_(j), n) * Leg->Val_(y.Val_(j), m) * _
Leg->Val_(x.Val_(j), p) * Leg->Val_(y.Val_(j), q)
NEXT : w.Set_(i, k, sum / j)
XYR_INDEX(p, q)
NEXT
sum = 0.0
FOR j = 0 TO az ' right side
sum += Leg->Val_(x.Val_(j), n) * Leg->Val_(y.Val_(j), m) * _
V->Val_(j, Cy)
NEXT : r.Val_(i) = sum / j
XYR_INDEX(n, m)
NEXT
Res = r / w ' solve system
IF StdDev THEN StdDev = Abw_(V, C1, C2, Cy)
END CONSTRUCTOR
DESTRUCTOR XYRegression()
IF Leg THEN DELETE Leg
END DESTRUCTOR
FUNCTION XYRegression.Val_(BYVAL X AS LA_S, BYVAL Y AS LA_S) AS LA_S
VAR xx = (X - Xmin) * Dx, yy = (Y - Ymin) * Dy, r = 0.0
VAR n = 0, m = 0
FOR i AS INTEGER = 0 TO Lines
r += Res.Val_(i) * Leg->Val_(xx, n) * Leg->Val_(yy, m)
XYR_INDEX(n, m)
NEXT : RETURN r
END FUNCTION
FUNCTION XYRegression.Abw_(BYVAL V AS LA_M PTR, _
BYVAL ChX1 AS UINTEGER = 0, _
BYVAL ChX2 AS UINTEGER = 1, _
BYVAL ChY AS UINTEGER = 2) AS LA_S
VAR i = V->Col_()
IF i < ChX1 ORELSE _
i < ChX2 ORELSE _
i < ChY THEN Err_ = "column missmatch" : RETURN 0.0
DIM AS LA_S sum = 0.0
FOR i = 0 TO V->Row_()
sum += (V->Val_(i, ChY) - Val_(V->Val_(i, ChX1), V->Val_(i, ChX2))) ^ 2
NEXT : RETURN sum / i
END FUNCTION
Anwendung
- Zur Anwendung von XYRegression.bas sind zwei Ergänzungen notwendig: Legendre Polynome und libFBla.
- Mit VAR test = XYRegression(Order, @V(), kx1, kx2, ky, 1, x1Min, x1Max, x2Min, x2Max) wird unter dem Namen test eine Fitting-Funktion generiert, welche den Meßwerten in V() möglichst nahe kommt:
- Der Parameter Order bestimmt die Ordnung der Ansatzfunktionen und damit die maximale Anzahl von Extrema im betrachteten Gebiet. Z. B. hat das Legendre Polynom der Ordnung 4 auch genau 4 Nullstellen und kann damit bis zu 3 Extrema abbilden. Die Ordnung 4 bei der Dimensionierung bedeutet also, dass maximal (4 - 1) ^ 2 = 9 Extrema auf der betrachteten Fläche abgebildet werden können. Je größer die Ordnung, desto besser können wellige Flächen angenähert werden. Aber desto größer ist auch die Rechenzeit und der Speicherbedarf (= 8 * (Ordnung + 1) ^ 2 nach erfolgter Dimensionierung). Bei überdimensionierter Ordnung nimmt die Ungenauigkeit sogar wieder zu (, aufgrund von numerischen Beschränkungen, welche sich zu größeren Fehlern aufsummieren.)
- Der Parameter V() ist vom Typ LA_V (= VeKtor, Array) und enthält drei Spalten (= 2. Dimension): x1, x2, y. Die Zeilen (= 1. Dimension) enthalten jeweils das Wertetrippel eines Meßpunktes. (Achtung, jede Zeile des Feldes geht als Meßwert in die Berechnung ein, also auch leere Zeilen (mit noch nicht definierten Werten).)
- Die Parameter kx1 und kx2 sind vom Typ UINTEGER und geben die Spaltennummern für die X-Werte an (Vorgabewerte: kx1 = 0, kx2 = 1).
- Der Parameter ky ist vom Typ UINTEGER und gibt die Spaltennummer für die Y-Werte an (Vorgabewert: ky = 2).
- Die 1 nach den Kanalnummern bewirkt, dass bei der Berechnung der Polynomparameter auch die mittlere Quadratische Abweichung bestimmt wird und mit test.StdDev ausgegeben werden kann. (Vorgabewert ist 0.)
- Die Parameter x1Min, x1Max und x2Min, x2Max bestimmen den Gültigkeitsbereich der Fitting-Funktion. Außerhalb dieses Bereiches liefert die Fitting-Funktion stark abweichende Werte. Der gewünschte Gültigkeitsbereich kann entweder vorgegeben werden (, wodurch z. B. auch Meßdaten extrapoliert werden können). Oder diese Parameter werden weggelassen und der Gültigkeitsbereich wird aus den Wertetrippeln des Feldes V() ermittelt.
- Die Funktion test.Abw_(V()) berechnet die mittlere quadratische Abweichung an den Meßpunkten x1, x2 zwischen den Werten der Fitting-Funktion und den Meßwert y. Kleine Rückgabewerte sind ein Zeichen für eine gute Näherung durch die Fitting-Funktion. Üblicherweise wird diese Funktion mit denselben Meßwerten V() aufgerufen, welche auch bei der Dimensionierung verwendet wurden.
- Ein Wert der Fitting-Funktion wird mit test.Val_(x1, x2) berechnet.
Beispiel
Beim Programmstart wird eine Fitting-Funktion vom Grad 4 berechnet. Nachdem das Bild aufgebaut wurde und die Meßpunkte sichtbar sind, kann mit den Tasten 0 bis 6 eine Neuberechnung mit entsprechendem Polynomgrad gestartet werden (Esc = Abbruch).
' This is file "XYRegressionTest.bas"
' (C) GPLv3 by Thomas[ dot ]Freiherr{ at }gmx{ dot }net
#INCLUDE ONCE "XYRegression.bas"
' grazy function to generate some test values
FUNCTION fn(BYVAL x1 AS DOUBLE, BYVAL x2 AS DOUBLE) AS DOUBLE
VAR x = (x1 - 379) / 319, y = (x2 - 112) / 199
RETURN 3 * EXP(-(x * x + y * y)) * SIN(x * y + 3)
END FUNCTION
' prepare graphics window
CONST ScW = 639, ScH = 399
SCREENRES ScW + 1, ScH + 1
' Prepare some test values, V() has 3 coloumns: x1, x2, y
RANDOMIZE TIMER
CONST Anzahl = 69 ' play with this number to see the differences
VAR V = LA_M(Anzahl + 1, 3), kx1 = 0, kx2 = 1, ky = 2
VAR mx = V.Val_(0, ky), mn = mx
FOR i AS INTEGER = 0 TO Anzahl
VAR x1 = RND * ScW
VAR x2 = RND * ScH
VAR yy = fn(x1, x2)
V.Set_(i, kx1, x1)
V.Set_(i, kx2, x2)
V.Set_(i, ky, yy)
IF yy > mx THEN mx = yy ' check extrema (for coloring)
IF yy < mn THEN mn = yy
NEXT
VAR f = mx - mn : IF f THEN f = 15 / f ' scale of colored area
VAR Order = 4, OrderMax = 6
DO
' Compute regression analysis
VAR test = XYRegression(Order, @V, kx1, kx2, ky, 1, 0.0, ScW, 0.0, ScH)
IF LEN(test.Err_) THEN _ ' on error
?"Fehler: "; test.Err_ : EXIT DO
WINDOWTITLE "Order " & Order & _
": (Grey colors valid, others are inaccurat)"
FOR x AS INTEGER = 0 TO ScW ' draw regression analysis
FOR y AS INTEGER = 0 TO ScH
VAR col = 16 + CUINT((test.Val_(x, y) - mn) * f)
PSET (x, y), col
NEXT
NEXT
FOR i AS INTEGER = 0 TO Anzahl ' draw measurement points
CIRCLE (V.Val_(i, kx1), V.Val_(i, kx2)), 3, 12, , , ,F
NEXT
WINDOWTITLE "Order " & Order & _
": Standard deviation: " & test.StdDev
' wait for key press to compute next step / warten auf Tastendruck
SLEEP
SELECT CASE INKEY
CASE CHR(27) : EXIT DO
CASE "0" : Order = 0
CASE "1" : Order = 1
CASE "2" : Order = 2
CASE "3" : Order = 3
CASE "4" : Order = 4
CASE "5" : Order = 5
CASE "6" : Order = 6
CASE CHR(255, 75), CHR(255, 80)
Order -= 1 : IF Order < 0 THEN Order = OrderMax
CASE ELSE
Order += 1 : IF Order > OrderMax THEN Order = 0
END SELECT
LOOP
English
See English forum.
| Zusätzliche Informationen und Funktionen | |||||||
|---|---|---|---|---|---|---|---|
|
|||||||
Bearbeiten
Versionen