This example is about a polynomial regression, a form of linear regression in which the relationship between the independent variable x and the dependent variable y is modeled as an nth order polynomial. For details see wikipedia -- Polynomial regression

Polynomial regression fits a nonlinear relationship between the value of x and the corresponding conditional mean of y, denoted y = F(x). (Even though it's called linear regression.) The code computes the parameters Ai of the function

so that the deviation between the polynomial and the measurement values gets minimal.F(x) =

A0

+ A1 x

+ A2 x^2

+ A3 x^3

+ ...

+ An x^n

This polynomial regression curve can be used to interpolate between measurement values or to extrapolate the measurment curve. Also the curve can be derived due to the known polynomial parameters.

This example consist of two parts

- 1) the UDT code defining a universal TYPE structure (for use in your project) and some

2) example code on how to use it.

This code contains the mathematical calculation of the polynomial parameters. For the linear algebra the library libFBla (is neccessary and) used. Save this code as 'PolReg.bas':

Code: Select all

```
' This is file "PolReg.bas"
' (C) LGPLv2 by Thomas[ dot ]Freiherr{ at }gmx{ dot }net
#INCLUDE ONCE "libFBla.bas"
TYPE Polynom
AS LA_S StDev
AS LA_V PolPara
DECLARE CONSTRUCTOR(BYVAL Order AS UINTEGER, BYVAL M AS LA_M PTR, _
BYVAL ChX AS UINTEGER = 0, BYVAL ChY AS UINTEGER = 1)
DECLARE PROPERTY Val_(BYVAL X AS LA_S) AS LA_S
END TYPE
' computes a polynom for data in matrix V
' ChX = column of X positions (default = 0)
' ChX = column of measurement values (default = 1)
' berechnet ein Polynom vom Grad Order aus den Werten in Matrix V
' ChX = Spalte der X-Positionen (Default = 0)
' ChY = Spalte der Meßwerte (Default = 1)
CONSTRUCTOR Polynom(BYVAL Order AS UINTEGER, BYVAL V AS LA_M PTR, _
BYVAL ChX AS UINTEGER = 0, BYVAL ChY AS UINTEGER = 1)
VAR az = V->Row_()
IF az < Order THEN ?"Error: not enough values" : EXIT CONSTRUCTOR
' prepare variables for the gradient equations
' coefficient matrix m and solution vector n
' Variablen fuer Berechnung der Normalengleichungen
VAR m = LA_M(Order + 1), n = LA_V(Order + 1)
' compute matrix (one triangle) and vector elements
' berechnen: Dreiecksmatrix und Lösungsvektor
FOR k AS INTEGER = 0 TO az
FOR i AS INTEGER = 0 TO Order
FOR j AS INTEGER = 0 TO i
m.Set_(j, i, m.Val_(j, i) + V->Val_(k, ChX) ^ (i + j))
NEXT
n.Val_(i) = n.Val_(i) + V->Val_(k, ChY) * V->Val_(k, ChX) ^ i
NEXT
NEXT
' mirror triangle matrix / Dreiecksmatrix spiegeln
FOR i AS INTEGER = 1 TO Order
FOR j AS INTEGER = 0 TO i - 1
m.Set_(i, j, m.Val_(j, i))
NEXT
NEXT
' compute the polynom parameters / Polynom-Parameter berechnen
PolPara = n / m
' compute the standard deviation / Standardabweichung berechnen
FOR i AS INTEGER = 0 TO az
StDev += (Val_(V->Val_(i, ChX)) - V->Val_(i, ChY)) ^ 2
NEXT
END CONSTRUCTOR
' comupte the value of the polynom at position X
' Berechnet den Wert des Polynoms an Position X
PROPERTY Polynom.Val_(BYVAL X AS LA_S) AS LA_S
WITH PolPara
VAR r = .Val_(0)
FOR i AS INTEGER = 1 TO .Ubound_()
r += .Val_(i) * X ^ i
NEXT : RETURN r
END WITH
END PROPERTY
```

- The polynomial parameters are calculated in the CONSTRUCTOR of the UDT Polynom. This gets executed by dimensioning ( VAR Pol = Polynom(Order, @m) ).
- The order of the polynomial is specified by the parameter Order.
- The matrix m of type LA_M (2D array) contains the input data. By default the X values are red from coloumn 0 and the Y values are red from coloumn 1. (Alternative coloumn numbers may be specified by parameters ChX (= channel X) and ChY, if required.
- The value of the polynomial at position x can be calculated by Pol.Val_(x).
- The polynomial parameters can be red by Pol.PolPara.Val_(n) with 0 <= n < Order.
- The standard deviation gets calculated in the CONSTRUCTOR and can be red by Pol.StDev.

Code: Select all

```
' This is file "PolReg_test.bas"
' (C) GPLv3 by Thomas[ dot ]Freiherr{ at }gmx{ dot }net
' Find more info at
' http://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)
' Weitere Informationen siehe
' http://de.wikipedia.org/wiki/Methode_der_kleinsten_Quadrate
#INCLUDE ONCE "vbcompat.bi"
' specify output format / Ausgabeformat spezifizieren
#DEFINE LA_FORMAT(_V_) FORMAT(_V_, "+#.0000e+000")
#DEFINE LA_TR !"\n"
' Try also SINGLE precision / Versuche SINGLE Berechnung
'TYPE AS SINGLE LA_S
'#DEFINE LA_Eps (1e-6)
#INCLUDE ONCE "PolReg.bas"
'#DEFINE english
' print out headers and values / Ausgabe Ueberschrift und Werte
#MACRO HEADER(_G_, _E_, _T_)
#IFDEF english
?:?_E_
#ELSE
?:?_G_
#ENDIF
?_T_
#ENDMACRO
' *** main ***
' here is the data input / Positionen X und Messwerte Y
VAR m = LA_M("1, 2" NL _
"2, 2.5" NL _
"3, 2.5" NL _
"4, 3.4" NL _
"5, 3.75" NL _
"6, 6.7" NL _
"7, 3")
' prepare output window / Ausgabefenster vorbereiten
CONST ScW = 639, ScH = 399
SCREENRES ScW + 1, ScH + 1, 8
WINDOWTITLE "PolRegTest.bas"
COLOR 0, 15
VAR Xm = 8, st = 2 * Xm / ScW
WINDOW (0, 0) - (Xm, 8)
VAR Order = 0 ' start order (simple average) / Start mit Mittelwert
DO
CLS
VAR Pol = Polynom(Order, @m)
HEADER("Polynomgrad:", _
"Polynomial order:" , Order)
HEADER("Polynom-Parameter:", _
"Polynomial parameters:", Pol.PolPara)
HEADER("Standardabweichung:", _
"Standard deviation:" , LA_FORMAT(Pol.StDev))
HEADER("Taste druecken!", _
"Press a key!", "")
' draw polynom / zeichne Polynomiale Ausgleichsfunktion
PSET (0, Pol.Val_(0)), 4
FOR x AS LA_S = st TO Xm STEP st
LINE - (x, Pol.Val_(x)), 12
NEXT
' draw measurement points / zeichne Messpunkte und Residuen
FOR i AS INTEGER = 0 TO m.Row_()
VAR x = m.Val_(i, 0), y = m.Val_(i, 1)
LINE (x, y) - (x, Pol.Val_(x)), 10 ' Residuen
CIRCLE (x, y), st, 9, , , , F ' Messpunkte
NEXT
' wait for key press to compute next step / warten auf Tastendruck
SLEEP
SELECT CASE INKEY
CASE CHR(27) : EXIT DO
CASE CHR(255, 75), CHR(255, 80)
Order -= 1 : IF Order < 0 THEN Order = m.Row_()
CASE ELSE
Order += 1 : IF Order > m.Row_() THEN Order = 0
END SELECT
LOOP UNTIL INKEY = CHR(27)
```

- The polynomial order must not be greater than the number of rows in the input data.
- Numerical issue may occur in case of high polynomial order or big number of input data (due to adding very big and very small values or due to big values exceeding the range of the variable type LA_S). To avoid this the input may be scaled to the interval [-1, 1] and/or the calculation of the gradient equations in the CONSTRUCTOR may be executed in several sub-sums.
- This regression can be used with 1D input only: y = F(x). For 2D input y = F(x, y) special (orthogonal) polynomials have to be used, see X / Y regression analysis.

http://www.freebasic-portal.de/code-bei ... n-231.html

Edit: link to libFBla added.