Code-Beispiel

X/Y-Regressionsanalyse

Lizenz:	Erster Autor:	Letzte Bearbeitung:
k. A.	TJF	08.11.2011

Anwendungsbeispiel Konturliniengrafik

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

Das Code-Beispiel wurde am 29.04.2011 von TJF angelegt.
Die aktuellste Version wurde am 08.11.2011 von TJF gespeichert.

Bearbeiten

Versionen