Polynomial regression curve (linear regression)

[TJF]

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

F(x) =
A0
+ A1 x
+ A2 x^2
+ A3 x^3
+ ...
+ An x^n

so that the deviation between the polynomial and the measurement values gets minimal.

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.

1) UDT code

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

Usage:

• 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.

2) Example usage

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)

Hints:

• 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.

Check also for updates:

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

Edit: link to libFBla added.

[dodicat]

Hi TJF

I've got this running, I had to google for your lib file.

I've got a folder on my desktop now called TJF to hold your includes.

I had planned to use another folder name, but my wife also uses this computer, so I reneged at the last minute.

It runs well, and I altered some y values to test.

Personally, I think that a freebasic array suffices to hold matrices/vectors, and that a matrix type just adds complexity.

Many years ago, on a small hand held programmable calculator, which could only handle an array of one dimension, I made a Gauss Jordan system solver.
I transferred it to BBC basic then Freebasic, and just for the hell of it I introduced complex elements (very easily due to Freebasic's operator capability).

It a little slower than the factorisation method, but it works ok.
I didn't bother to extract the determinant, for I didn't really need it.

Full pivoting performed at every stage, for rounding errors were a problem on the little calculator.

Code: Select all

'GAUSS JORDAN FOR REAL OR COMPLEX SOLUTIONS
'define the variables
type complex
    as double re,im
end type
operator +(n1 as complex,n2 as complex) as complex
dim as complex n
n.re=n1.re+n2.re
n.im=n1.im+n2.im
return n
end operator

operator -(n1 as complex,n2 as complex) as complex
dim as complex n
n.re=n1.re-n2.re
n.im=n1.im-n2.im
return n
end operator

operator *(n1 as complex,n2 as complex) as complex
dim as complex n
n.re=n1.re*n2.re - n1.im*n2.im
n.im=n1.im*n2.re + n1.re*n2.im
return n
end operator

operator /(n1 as complex,n2 as complex) as complex
dim as complex n
dim as double d
d=n2.re^2+n2.im^2
n.re=(n1.re*n2.re+n1.im*n2.im)/d
n.im=(n1.im*n2.re - n1.re*n2.im)/d
return n
end operator

'GAUSS-JORDAN METHOD FOR REAL OR COMPLEX LINEAR EQUATIONS
Sub gauss(mat() As complex,rhs() As complex,solution() As complex)
Dim dd As Long:dd=Ubound(rhs)
Redim _a(1 To (dd-1 + 2) * (dd-1 + 1)) As complex 'main working array 
Redim mat2(1 To dd,1 To dd+1) As complex 'holds matrix plus rhs
dim null as complex
For rows As Long=1 To dd
    For cols As Long=1 To dd
        mat2(rows,cols)=mat(rows,cols)  'TRANSFER to a copy
    Next cols
Next rows

'add the rhs to the copy
For fill As Long =1 To dd
    mat2(fill,dd+1)=rhs(fill)
Next fill

Dim count As Long
 For rows As Long=1 To dd
  For cols As Long=1 To dd+1   
count=count+1
_a(count)=mat2(rows,cols) 'Fill the working array
 Next cols
Next rows

    Dim d As Long:d=dd-1
Dim As Long x,v,z,w,dimcount,vvcount,k2,vv,x1,k1
Dim y As complex
 v = 1: w = d + 1: x = 0: dimcount = 0
 
'   PIVOT THE MATRIX BY REAL 
For x1=dimcount To dd-1
    For k1=x1+1 To dd-1
      If Abs(_a(1+k1*(dd+1)+dimcount).re) > Abs(_a(1+x1*(dd+1)+dimcount).re) Then
For g As Long=0 To dd
    swap _a(1+k1*(dd+1)+g),_a(1+x1*(dd+1)+g)
    Next g
          End If 
       Next k1
   Next x1:dimcount=dimcount+1 
'  END PIVOT 

Redim mult(1 To w + 1) As complex
vvcount = 0: k2 = 0
twenty:If vvcount = 0 Then
y=_a(v)
 For z = v To v + w
    _a(z) = _a(z) / y 
   Next z 'across each row  
End If
If vvcount > 0 Then 
   mult(vvcount) =null -_a(vv) 
end if

 If w = 1 Then Goto eighty
 v = v + d + 2
If v > ((d+1)*(d+2)) Then:v=(d+3)+(d+3)*(d+1-w):Goto fifty:End If
vv = v: vvcount = vvcount + 1
 Goto twenty                        'all rows
fifty: vvcount = 1 + k2
For z = v To v + w
 _a(z) = _a(z) + mult(vvcount) * _a(1 + (d + 3) * (d + 1 - w) + x)
 x=x+1
 Next z
 v = v + d + 2: x = 0
 If v > ((d + 1) * (d + 2)) Then
w = w - 1
Redim mult(1 To w + 1) As complex
v = (1 + (d + 3) * (d + 1 - w))

If w<>1 Then '  FULL RE-PIVOT AT EACH STAGE BY REAL ,(inline for speed) 
For x1=dimcount To dd-1
    For k1=x1+1 To dd-1
      If Abs(_a(1+k1*(dd+1)+dimcount).re) > Abs(_a(1+x1*(dd+1)+dimcount).re) Then
For g As Long=0 To dd
    Swap _a(1+k1*(dd+1)+g),_a(1+x1*(dd+1)+g)
    Next g
          End If 
         Next k1
     Next x1:dimcount=dimcount+1 
Endif
'END RE-PIVOT
vvcount = 0: k2 = 0:Goto twenty                           'all rows
End If
 k2 = k2 + 1: Goto fifty
eighty:   Redim  solution(1 To d + 1) As complex 'make the answer array
For w = 1 To d + 1
solution(w) = _a(1 + (d + 3) * (d + 1 - w) + w)  'back substitute

For x = 2 To w
    solution(w) = solution(w) - (_a(1 + (d + 3) * (d + 1 - w) + w - (x - 1)) * solution(x - 1))
Next x
Next w
'reverse the array to get the answer
Dim As Long lb,ub:lb=Lbound(solution):ub=Ubound(solution)
    For n As Long=Lb To Int((lb+Ub)/2):Swap solution(n),solution(ub+lb-n):Next
End Sub

' NOT PART OF GAUSS-JORDAN, used to check the answer
Sub mv(m1() As complex,m2() As complex,ans() As complex) 'MATRIX x VECTOR
   Dim s As complex
   dim n as double=ubound(m1)
   redim ans(1 to n)
    For i As Integer=1 To n
        s.re=0:s.im=0
        For k As Integer = 1 To n
            s=s+m1(i,k)*m2(k)
        Next k
        ans(i)=s
        Next i
    End Sub
' ***************** EXAMPLE ****************************    
redim as complex mat(0,0),rhs(0)

dim di as integer 
''print "FOR COMPLEX ENTRIES, REAL a comma then IMAGINARY"
print "EG.  -9,12 "
print
Input "enter dimension for n x n equations i.e. what is n   ",di
if di<=0 then end

Redim mat(1 To di,1 To di) 'make a matrix (1 to di,1 to di) to fill
Redim rhs(1 To di)        'make a rhs vector (1 to di) to fill
'PLEASE NOTE, ANSWER VECTORS MUST BE REDIM
Redim  answer(0) as complex 'make two empty arrays for the answers

redim answer2(0) as complex 'this one for the check multiply
'INPUT THE MATRIX
For rows As Integer=1 To di
    For cols As Integer=1 To di
         ''Print " row ";rows;" column ";cols;" =      "; 'fill it up
         ''Input mat(rows,cols).re,mat(rows,cols).im
         mat(rows,cols)=type<complex>(rnd*1-rnd*1,rnd*2)
    Next cols
   '' print
Next rows
print "Random matrix made"
print
' INPUT THE RIGHT HAND SIDE VECTOR
print "NOW THE RIGHT HAND SIDE VECTOR"
For rows As Integer=1 To di
    ''Print " RHS row ";rows;" =      ";      'fill it up
    ''Input rhs(rows).re,rhs(rows).im
    rhs(rows)=type<complex>(rnd*5-rnd*5,rnd*4)
    print rhs(rows).re,rhs(rows).im;" *i "
Next rows

gauss(mat(),rhs(),answer())  'GET answer()

print
print "SOLUTION"
for z as integer=1 to ubound(answer)
    print answer(z).re,answer(z).im;" *i "
next z
print 

mv(mat(),answer(),answer2())   'RE MULTIPLY TO CHECK

print "CHECK:"
print "MULTIPLY THE SOLUTION BY THE MATRIX GIVES BACK THE RIGHT HAND SIDE"
print
for z as integer=1 to ubound(answer2)
    print answer2(z).re,answer2(z).im;" *i "
next z
print "DONE"
sleep

[TJF]

Hello dodicat!

I've got this running, I had to google for your lib file.

Sorry for your trouble and thanks for the hint! I fixed it by adding a link in the OP.

I've got a folder on my desktop now called TJF to hold your includes.

Glad to read that. I'll do my best to fill the folder :)

Personally, I think that a freebasic array suffices to hold matrices/vectors, and that a matrix type just adds complexity.

Yes, it adds complexity. But when you're used to it, the source code gets clearer. And your coding gets faster because you need not care about the index. To copy a complete matrix it's just VAR b = a. No DIM, no FOR : ... : NEXT, etc.

... just for the hell of it I introduced complex elements ...

You can use your complex UDT with libFBla. Just start your code with (for libFBla version >= 0.2)

• TYPE LA_S AS complex
  #INCLUDE ONCE "libFBla.bas"

First you'll have to complete your complex type. Add CONSTRUCTORS and some more OPERATORs (like CAST, +=, ...) and overload some functions like SQR(). But when your complex UDT is a full featured number type you can base all libFBla functions on this type.

[Lost Zergling]

Hello,
dodicat wrote :

Personally, I think that a freebasic array suffices to hold matrices/vectors, and that a matrix type just adds complexity.

I agree.
Nevertheless, I see two weakpoints in FB (from a neophyte's point of view)(in maths programming) :
1 - The lack of an"easy" maths/vectors dedicated instruction set in FB (vs other languages) (vector/maths oriented i/s)
2 - The lack of automated memory management for arrays processing tasks (can lead to massive redim preserve sometimes inside types)
Many years ago (before lzle), I tried to program a Jordan Gauss pivot using strings lists, the goal was to handlle very long factorial decomposition on big numbers as strings on small matrices so as to try to "smart" simplificate fractions ("half symbolic" process)(ie : "3*7*5/13*11", and so on).
As a result, memory consumption and processing times have exploded because of the cases of large numbers which were nevertheless very close but without easily smartly detectable simplifications of fractions.
I hate maths,... ;-)

[ardmin]

Hi,

When I try to use polynomial regression with big numbers for x and double check the result in excel there are differences and seems that the program here is wrong. For instance keep the y data same and set x=1000, 1001, 1002, etc - then the error is in the 7th number. If you set x=30001, 30002, 30003 etc then the polynomial coefficients are absolutely different. What could be the reason for that and how to deal with it?

Thank you very much.

[TJF]

Hi ardmin!

It's hard to give good advice without knowing your source code.

... then the error is in the 7th number.

Computing in single precision is 6.9 digits accuracy (23 bits mantissa, see IEEE 754). Double precision has

• a bigger range (11 bits exponent)
• more precision (52 bits mantissa -> 15.9 digits)

Perhaps you can solve both issues by computing in double precision. Start your code by

Code: Select all

TYPE AS DOUBLE LA_S
#INCLUDE ONCE "libFBla.bas"

Regards

[ardmin]

Hi TJF,

Thank you for your reply. I used the source code above, just put values x=10001, 10002, etc. y=1.001, 1,002, etc. and used degree=2.
Next I compared the results with excel and it doesn't match at all.

Best,
ardmin

[TJF]

Hi ardmin!

It's hard to give good advice without knowing your source code.

Also mind the hints in the first post!

Regards