Polynomial curve fitting

[mrcoollion]

I am trying to build the best polynomial curve fitting lines/formula's from a number of X Y points
Has anybody done this before with LC ? Are there mathematical functions in LiveCode I can use for this (like in excel with LINEST or POLIFIT in Mathlab).
Or does anybody has some sample script available?

Regards,

Paul

[richmond62]

Presumably you would have a polygon graphic
and you would set the points of your polynomial curve to that.

A few years ago I made a "spirograph" thing in LiveCode that did just that:
can't for the life of me think where it might be at the moment.

[mrcoollion]

Presumably you would have a polygon graphic
and you would set the points of your polynomial curve to that.

A few years ago I made a "spirograph" thing in LiveCode that did just that:
can't for the life of me think where it might be at the moment.

Thanks for the reply R...
However, in my case I have many points (X,Y coordinates) and I want to get the best fitting line through those points.
This means I want to create a polynomial function that approximates a set of data points

graphpoints.png (5.19 KiB) Viewed 9150 times

PointsGraph2.gif (10.6 KiB) Viewed 9147 times

[richmond62]

Presumably the polynomial function will "spit out" coordinates for your curve which can then
be used as the points for a graphic object.

[mrcoollion]

Presumably the polynomial function will "spit out" coordinates for your curve which can then
be used as the points for a graphic object.

More or less...
I will be using the polynomial function to predict the next coordinate.

[richmond62]

Pretty crude, but you get the point . . . .

Erm . . . hang on a mo' . . . the points.
-

[dunbarx]

Hermann.

Get in here.

Craig

[[-hh]]

First please read https://en.wikipedia.org/wiki/Polynomial_regression
to understand that polynomial regression is a special case of multiple linear regression.

So all you essentially need is a fast handler for computing the inverse of the matrix X'X where X is given in the wiki above.

This is basic linear algebra and not hard to implement in LC.
[See for example the last code snippet here: viewtopic.php?p=181709#p181709]

Note.
When computing polynomial fits for a degree > 3 you might run into rounding errors with LC. You could then use a bignum library of javascript via a browser widget or of python via shell).

[mrcoollion]

I almost wracked my brain on it and it took me 2 days, but I made a working routine that calculates the Polynomial Coefficients up to a maximum of 7 degrees.
Is anybody interested in it also for testing?
I do think that LC should have some more mathematical functions available like this one, as most development languages have. Building this was not good for my blood pressure...

Let me know and I will post the code.

[[-hh]]

Why do you ask if you already have a working routine? What a waste of time ...

[mrcoollion]

I did not have a working routine.. I just build it.

[dunbarx]

Hermann.

MrCoolion was in the process of working on his routine when he asked if anyone had already done something similar. This is not unusual. I have often stopped in the middle of a task to see if there is an existing solution somewhere.

Craig

[mrcoollion]

Thanks Craig,

You gave me back the positive vibe of today. I will post my code in a few moments for anyone to use and audit .

Regards,

Paul

[mrcoollion]

Here is the command with which one can determine the Polynomial Coefficients of a formula for curve fitting.
Copy>Paste the complete code in one button and have a look at the array OutPolCo at the breakpoint.

Command consists of following routines.

• Build a polynomial regression matrix for Xi (tXMatrix)

• Build the XiYi matrix (tYMatrix)

• Combine the Xi and XiYi matrixes into an augmented matrix for Gauss-Jordan elimination (tXMatrix)

• Do the Gauss-Jordan elimination on the augmented matrix (result is in taSolution)

Have fun with it..... And if you have any improvement please let me know and post them here.

Code: Select all

on mouseup
   // See for conformation and information: https://arachnoid.com/sage/polynomial.html
   put "-1,0,1,2,3,5,7,9" into InXCoordinates
   put "-1,3,2.5,5,4,2,5,4" into InYCoordinates
   put 4 into InDegrees
    ---- More test data
   --put "0,1,2,3,4,5,6,7,8,9,10" into InXCoordinates  // 4 test
   --put "1,6,17,34,57,86,121,162,209,262,321" into InYCoordinates  // 4 test
   --put 2 into InDegrees  / 4 test
   ------
   CalcPolynomialCoefficients InDegrees,InXCoordinates,InYCoordinates,OutPolCo
   ------
   breakpoint // See OutPolCo for the coefficients
end mouseup

command CalcPolynomialCoefficients InDegrees,InXCoordinates,InYCoordinates @OutPolCo
   // See for conformation and information: https://arachnoid.com/sage/polynomial.html
   --put "-1,0,1,2,3,5,7,9" into InXCoordinates // 4 test
   --put "-1,3,2.5,5,4,2,5,4" into InYCoordinates // 4 test
   -- put 4 into InDegrees // test
   // Higest Aray number in OutPolCo gives the coefficient for the higest degree of X in a equation like y= a*X4+b*X3+c*X2+d*X+e
   // In this case with the test numbers OutPolCo[5] is the coefficient for X4
   -----------------------------------------------------------------------------------
   put the number of items of InXCoordinates into tNumberOfXItems
   put the number of items of InYCoordinates into tNumberOfYItems
   if tNumberOfXItems <> tNumberOfYItems
   then
      answer "Amount of X and Y coordinates is not equal!"
      exit CalcPolynomialCoefficients
   end if
   if InDegrees > 7
   then
      answer "Max number of degrees is 7. Larger polynomial degrees risk exceeding the resolution of the variables in use !"
      exit CalcPolynomialCoefficients
   end if
   if tNumberOfXItems < InDegrees+1
   then
      answer "Need a minimum of  " & InDegrees+1 & " sets of x,y coordinates!"
      exit CalcPolynomialCoefficients
   end if
   ---- fill X matrix -----------------------------------------
   repeat with MVC = 1  to  InDegrees+1
      Repeat with MHC = 1  to  InDegrees+1
         if MVC=1 and MHC=1
         then
            put tNumberOfXItems into tSumXiP 
         else
            put MVC+MHC-2 into tPower
            repeat for each item tXi in InXCoordinates
               --
               put "put tXi^"&tPower&" into tXiP" into tScriptLine
               do tScriptLine
               --
               add tXiP to tSumXiP
            end repeat
         end if
         put tSumXiP into tXMatrix[MVC,MHC] 
         put 0 into tSumXiP
      end Repeat
   end repeat
   -- Fill Y Matrix (MVC stands for Matrix Vertical Column)----
   repeat with MVC = 1  to  InDegrees+1
      put 0 into tCounter
      repeat for each item tYi in InYCoordinates
         add 1 to tCounter
         if MVC=1
         then
            add tYi to tSumXiPYi
         else
            put MVC-1 into tPower
            put item tCounter of InXCoordinates into tXi
            --
            put "put tXi^"&tPower&" into tXiP" into tScriptLine
            do tScriptLine
            --
            add tXiP*tYi to tSumXiPYi
         end if
      end repeat
      put tSumXiPYi into tYMatrix[MVC]
      put 0 into tSumXiPYi
   end repeat
   -- Combine matrixes by adding Y matix as last column to X matrix 
   -- creates an augmented matrix
   repeat with MVC = 1  to  InDegrees+1
      put tYMatrix[MVC] into tXMatrix[MVC,InDegrees+2]
   end repeat
   ----------------------------------------
   // Gauss-Jordan elimination 
   put InDegrees+1 into tNbrRows 
   put InDegrees+2 into tNbrColumns 
   --
   repeat with tColumnNbr = 1  to  tNbrRows 
      repeat with tRowNbr = 1  to  tNbrRows 
         if tColumnNbr <> tRowNbr
         then
            put tXMatrix[tRowNbr,tColumnNbr]/tXMatrix[tColumnNbr,tColumnNbr] into tCNumber
            repeat with k =1 to tNbrColumns
               put tXMatrix[tRowNbr,k]-tCNumber*tXMatrix[tColumnNbr,k] into tXMatrix[tRowNbr,k]
            end repeat
         end if
      end repeat
   end repeat
   --
   repeat with tRowNbr = 1  to  tNbrRows
      put tXMatrix[tRowNbr,tNbrRows+1]/tXMatrix[tRowNbr,tRowNbr]into taSolution[tRowNbr]
   end repeat
   ------------------------------------------
   put taSolution into OutPolCo
end CalcPolynomialCoefficients

[[-hh]]

This is the well known polysolve algorithm (which is the general linear model solution specialized for the polygonial hypothesis).

See an online computation (incl. your example and a graphical output) for example here:
https://arachnoid.com/polysolve/index.html

Now you translated that algorithm to LC, fine, well done!

Sadly, the problem with that (and generally with using LC for advanced statistical computations) is that you may get tremendous rounding errors when using floating point math (as LC does), depending on the dataset and the degree of the polynomial.
This is for example described in "popular" math wording on the above page at the bottom in section "Technical notes and limits".

That's why some people who need that don't use LC but a program that can handle bignums (as mathematica, matlab, sage, ..., python).