Polynomials [Resolved]

[yrwyddfa]

I need to be able to calculate a polynomial trendline. I've cracked Gauss-Jordan elimination (and written the code) to solve the equations - but how do I generate the equations in order to solve them.

I'm a bit simple so some explanation will be required.

. . . I know there's a genius out there somewhere . . . .

[wossname]

You could use a genetic algorithm. It's an interesting challenge but very fast for high order polynomials.

[yrwyddfa]

[deleted because it was obsolete]

[yrwyddfa]

Here's my effort: it appears to work identically to Excel:

Option Explicit

Public Type POINT
    x As Double
    y As Double
End Type


Public Function Trend(Data() As POINT, ByVal Degree As Long) As POINT()

    'degree 1 = straight line y=a+bx
    'degree n = polynomials!!
    
    Dim a() As Double
    Dim Ai() As Double
    Dim B() As Double
    Dim P() As Double
    Dim SigmaA() As Double
    Dim SigmaP() As Double
    Dim PointCount As Long
    Dim MaxTerm As Long
    Dim m As Long, n As Long
    Dim i As Long, j As Long
    Dim Ret() As POINT
    
    Degree = Degree + 1
    
    MaxTerm = (2 * (Degree - 1))
    PointCount = UBound(Data) + 1
    
    ReDim SigmaA(MaxTerm - 1)
    ReDim SigmaP(MaxTerm - 1)
    
    ' Get the coefficients lists for matrices A, and P
    For m = 0 To (MaxTerm - 1)
        For n = 0 To (PointCount - 1)
            SigmaA(m) = SigmaA(m) + (Data(n).x ^ (m + 1))
            SigmaP(m) = SigmaP(m) + ((Data(n).x ^ m) * Data(n).y)
        Next
    Next
    
    ' Create Matrix A, and fill in the coefficients
    ReDim a(Degree - 1, Degree - 1)
    For i = 0 To (Degree - 1)
        For j = 0 To (Degree - 1)
            If i = 0 And j = 0 Then
                a(i, j) = PointCount
            Else
               a(i, j) = SigmaA((i + j) - 1)
            End If
        Next
    Next
    
    ' Create Matrix P, and fill in the coefficients
    ReDim P(Degree - 1, 0)
    For i = 0 To (Degree - 1)
        P(i, 0) = SigmaP(i)
    Next
    
    ' We have A, and P of AB=P, so we can solve B because B=AiP
    Ai = MxInverse(a)
    B = MxMultiplyCV(Ai, P)
    
    ' Now we solve the equations and generate the list of points
    PointCount = PointCount - 1
    ReDim Ret(PointCount)
    
    ' Work out non exponential first term
    For i = 0 To PointCount
        Ret(i).x = Data(i).x
        Ret(i).y = B(0, 0)
    Next
    
    ' Work out other exponential terms including exp 1
    For i = 0 To PointCount
        For j = 1 To Degree - 1
            Ret(i).y = Ret(i).y + (B(j, 0) * Ret(i).x ^ j)
        Next
    Next
    
    Trend = Ret
    
End Function

Public Function MxMultiplyCV(Matrix1() As Double, ColumnVector() As Double) As Double()

    Dim i As Long
    Dim j As Long
    Dim Rows As Long
    Dim Cols As Long
    Dim Ret() As Double
    
    Rows = UBound(Matrix1, 1)
    Cols = UBound(Matrix1, 2)
    
    ReDim Ret(UBound(ColumnVector, 1), 0) 'returns a column vector
    
    For i = 0 To Rows
        For j = 0 To Cols
            Ret(i, 0) = Ret(i, 0) + (Matrix1(i, j) * ColumnVector(j, 0))
        Next
    Next
    
    MxMultiplyCV = Ret
    
End Function

Public Function MxInverse(Matrix() As Double) As Double()
    
    Dim i As Long
    Dim j As Long
    Dim Rows As Long
    Dim Cols As Long
    Dim Tmp() As Double
    Dim Ret() As Double
    Dim Degree As Long
    
    Tmp = Matrix
    
    Rows = UBound(Tmp, 1)
    Cols = UBound(Tmp, 2)
    Degree = Cols + 1
    
    'Augment Identity matrix onto matrix M to get [M|I]
    ReDim Preserve Tmp(Rows, (Degree * 2) - 1)
    For i = Degree To (Degree * 2) - 1
        Tmp(i Mod Degree, i) = 1
    Next
    
    ' Now find the inverse using Gauss-Jordan Elimination which should get us [I|A-1]
    MxGaussJordan Tmp
    
    ' Copy the inverse (A-1) part to array to return
    ReDim Ret(Rows, Cols)
    For i = 0 To Rows
        For j = Degree To (Degree * 2) - 1
            Ret(i, j - Degree) = Tmp(i, j)
        Next
    Next
    
    MxInverse = Ret
    
End Function

Public Sub MxGaussJordan(Matrix() As Double)
    
    Dim Rows As Long
    Dim Cols As Long
    Dim P As Long
    Dim i As Long
    Dim j As Long
    Dim m As Double
    Dim d As Double
    Dim Pivot As Double
    
    Rows = UBound(Matrix, 1)
    Cols = UBound(Matrix, 2)

    ' Reduce so we get the leading diagonal
    For P = 0 To Rows
        Pivot = Matrix(P, P)
        For i = 0 To Rows
            If Not P = i Then
                m = Matrix(i, P) / Pivot
                For j = 0 To Cols
                    Matrix(i, j) = Matrix(i, j) + (Matrix(P, j) * -m)
                Next
            End If
        Next
    Next
    
    'Divide through to get the identity matrix
    'Note: the identity matrix may have very small values (close to zero)
    'because of the way floating points are stored.
    For i = 0 To Rows
        d = Matrix(i, i)
        For j = 0 To Cols
            Matrix(i, j) = Matrix(i, j) / d
        Next
    Next
    
End Sub

Enjoy!

[WonkoTheSane]

I was imagining you wanted to construct the actual equation that represents the trend line.

Something like "y = 3.2x^4 + 2.923x^3 + 0.012x^2 - 71.322x + 4".

[yrwyddfa]

You can derive the relevant expression from the code posted. Look at the last two loops of the Trend function.

[wossname]

Originally posted by yrwyddfa
You can derive the relevant expression from the code posted. Look at the last two loops of the Trend function.

Do you intend to add this capability? It would be very useful to me if you would be willing to share it.

[yrwyddfa]

    Equation = "y=" & Format$(B(0, 0), "0.00000") & " + "
    For j = 1 To Degree - 1
        Equation = Equation & Format$(B(j, 0), "0.00000") & "x^" & j & " + "
    Next
    Equation = Left$(Equation, Len(Equation) - 3)
    
    Debug.Print Equation

Add this to the bottom of the Trend function.

[wossname]

Cool.

Could it generate this red curve (or an approximation) from the noisy data?

[yrwyddfa]

Yup; that looks like quite a high order polynomial.

What I've done (and I'm not sure the definition is correct) but the degree parameter semantics are:

1:y=a+bx
2:y=a+bx+bx^2
3:y=a+bx+bx^2+bx^3

I imagine to generate that curve you would need a 4, or 5.

If you don't mind me asking . . .. what field data did you get a sigmoidal curve from?

[wossname]

I was taught

y=ax^5 + bx^4 + cx^3 + dx^2 + ex + f


But anyway, I generated that data at random using a simple y=x^3+(+/-RandomNumber).

Then I warped it a bit in a photo editor.

How would I call your Trend() function to process this data?

Lets say Raw(0-999) contains the point data.

[wossname]

Originally posted by yrwyddfa
Yup; that looks like quite a high order polynomial.

I would'nt call 4 or 5 high order though 20 yes. 30 certainly.

[yrwyddfa]

How would I call your Trend() function to process this data?

Lets say Raw(0-999) contains the point data.

Dim t() As POINT
    Dim P(33) As POINT
 
 
    P(0).y = 2
    P(1).y = 4
    P(2).y = 3
    P(3).y = 6
    P(4).y = 5
    P(5).y = 7
    P(6).y = 8
    P(7).y = 6
    P(8).y = 9
    P(9).y = 8
    P(10).y = 7
    P(11).y = 6
    P(12).y = 5
    P(13).y = 4
    P(14).y = 3
    P(15).y = 2
    P(16).y = 3
    P(17).y = 4
    P(18).y = 5
    P(19).y = 4
    P(20).y = 5
    P(21).y = 6
    P(22).y = 7
    P(23).y = 8
    P(24).y = 90
    P(25).y = 9
    P(26).y = 8
    P(27).y = 7
    P(28).y = 8
    P(29).y = 98
    P(30).y = 7
    P(31).y = 67
    P(32).y = 56
    P(33).y = 45
 
    t = Trend(P, 5)

Also don't forget to fill in the 'x' axis (poss with a for loop?)

The array t will then hold the x, and y coordinates for the line which you can then plot.

Let me know how you get on . . .

[yrwyddfa]

y=ax^5 + bx^4 + cx^3 + dx^2 + ex + f

Well, it's identical.

the equation

of a line is: y = a+bx
of a parabola is y=a+bx+bx^2

etc etc.

Your representation is simply back to front from mine; they both equate to the same thing,

[yrwyddfa]

Is this going to end being sanded?

[wossname]

Not in the near future, no. I have another lib in the making (nearing completion) and this would fit brilliantly.

I already have it graphically plotting an equation given the "Y=blahablablah" which is pretty fast (compiles it internally before running repeated evaluations) and this would be cool to include.

I won't need the Trend function to return a point() array, just the string y=...

Then I can plug that string back into the plotter and it should draw the trend alongside the datasample. Cool beans.



I'll try this out tonight. Cheers

[yrwyddfa]

Good stuff Let me know how you get on.

[wossname]

OMG!!

Dude this RULES. I can scarcely contain myself . I had to take a snap for posterity...



The Graphic calulator is showing the "REAL" best fit curve (that is, the curve that the random data was created from). So it's safe to say it is the best fit for this data.

The laptop is my program with your code added. The squiggly black lines are lines between data points, I don't know how to plot points without lines yet.

The red curve was generated by your algorithm based on the squiggly line data points. Its as close to dead on perfect as I could have hoped for.

yrwyddfa, congrats of a whupass algorithm mate! Top quality. I'm going to optimise it for .net and see how fast I can get it running.

[wossname]

ph34r yrwyddfa!

The attached image is a plot screenshot of

• the raw data points (squiggly black lines)
• the Trend() in RED
• the "Real" best fit curve (BLUE) from which the random raw data was created.

You will have to zoom in to be able to see the red and blue clearly. That is because they overlap almost perfectly. THAT is cool.




Scary thing is that I only gave trend() 25 points to work with! I imagine the accuracy is infallible with 100 or more!

Good work Sir.

[yrwyddfa]

Pleased you liked it . . .

The algorithm itself is quite optimised because I used Gauss-Jordan reduction to get the inverse matrix.

If you use the classic multiplying method you would need n! calculations, but using Gauss-Jordan you need n^3.

If you are using polynomials of degree 4 or higher then this algorithm is efficient (in terms of number of calculations) So it IS ineffecient for straight lines - particularly there is a shortcut to get the inverse of 2x2 matrix.

If you are converting this to .Net you'll probably need to create a Matrix class with the following methods: Prop Get/Let Element, Inverse, GaussJordan, Multiply, Augment. This covers most of the major mathematics in the algorithm.

The worst part of it all is that I couldn't find ANY source code examples (VB or otherwise) that do this ANYWHERE on the web.

Maybe my googling is not up to it . . .

(FYI - the technique is called least-squares curve fitting)

[yrwyddfa]

What does ph34r mean?

[wossname]

Sorted out the dot drawing problem and optimised your code for .net...



I don't know anything about matrices (except the kind you sleep on ) though. Can you point me to any URLs that will tell me the theory, so i can use it to build a matrix class?

I'll give you a credit in the documentation for your contribution to the project, PM me if you want me to put your real name or any other info in also

This is really cool

[dglienna]

ph34r is l33t-speak for FEAR!

if that's your EEG, then you are one sick kitty!

[wossname]

yeah, but it looks sweet though right?

[yrwyddfa]

Originally posted by wossname
Sorted out the dot drawing problem and optimised your code for .net...



I don't know anything about matrices (except the kind you sleep on ) though. Can you point me to any URLs that will tell me the theory, so i can use it to build a matrix class?

I'll give you a credit in the documentation for your contribution to the project, PM me if you want me to put your real name or any other info in also

This is really cool

Any good book on discrete mathematics will do the job for matrices: the one I used is called 'A First Course in Discrete Mathematics' By Molluzo, and Buckley ISBN: 0-534-05310-6.

I have no idea if it's still in print. Trust me - you need the book as a constant reference when you start doing this stuff.

I never really found any excellent internet references on any of the stuff I used in this algorithm.

[yrwyddfa]

I'll give you a credit in the documentation for your contribution to the project, PM me if you want me to put your real name or any other info in also

Nah. Yrwyddfa@vbforums will do. Perhaps there is some way the mods can set up an area with v. useful algorithms so you could stick a URL in there as well.

Places such as these are useful for finding out how to do new stuff; I wouldn't want other people to go through days of research in order to find out how to stuff that is essentially public domain anyway.

[opus]

Cool code, but why do I get a different formula back if i run the function a second time?

[wossname]

Originally posted by opus
Cool code, but why do I get a different formula back if i run the function a second time?

Hmm, odd. Could you please post your test code?

[opus]

Public Sub test()
Dim data() As POINT
Dim trendwerte() As POINT
Dim i As Integer
ReDim data(2000)
ReDim trendwerte(2000)
'Read data
For i = 0 To 1999
    'Column 1 has just to X Values (for example 0,1,2........)
    data(i).x = ActiveWorkbook.Sheets(1).Cells(i + 2, 1).Value
    'Column 2 has the Y Values (like your example X^3 +/- randomvalue)
    data(i).y = ActiveWorkbook.Sheets(1).Cells(i + 2, 2).Value
Next i
trendwerte = Trend(data, 3)
For i = 0 To 1999
    'Column 3 recieves the trenddata
    ActiveWorkbook.Sheets(1).Cells(i + 2, 3).Value = trendwerte(i).y
Next i
 
End Sub

I'm using the code as posted (including the Debug.Print for the formula).
in 3 consecutive runs of the test Sub I get these values for the formula:
y=13720761,62709 + -45666,97498x^1 + -10,67540x^2 + 1,09310x^3
y=-16117336,06632 + -5163,26485x^1 + 37,59203x^2 + 1,02608x^3
y=9902854,69155 + 47807,57558x^1 + 2,81787x^2 + 0,87641x^3

[opus]

Sorry,
I found my error.
I used EXCEL to produce my data (including a call for Randomnumbers). and I used the trend-function in that file. Each time I run the trend Function, the value that are calculated by the RandomNumber changed. And for that reason the resulting formula has different values.

Me stupid

[dglienna]

Add the Checkmark and the word [RESOLVED] to the subject of the first post in your thread
if your question has been answered satisfactorily.

[wossname]

Originally posted by dglienna
Add the Checkmark and the word [RESOLVED] to the subject of the first post in your thread
if your question has been answered satisfactorily.

I bet you were the milk monitor at school weren't you?

[dglienna]

it looked like it needed *something*. next time, I'll just send a PM to have it manually [RESOLVED] with the tick.

[wossname]

.Net Framework graphics suck.

I'm going to need to redesign the drawing routines.

[yrwyddfa]

DC, and GDI then . . .

[wossname]

Might be OK after all. Had a big problem with scaling using the built-in .net matrices. I've got to scale all the coords manually from now on and use doubles until the last minute conversion to singles.

[azteched]

If you want the inverse of a matrix, there's faster ways than Gauss-Jordan - read it up in a computational maths textbook

[wossname]

Maths books make me ill. All those pointless arbitrary symbols. These math boffins claim to be so smart but they can't write stuff down so dumx0rs like me can understand it.

[azteched]

It rather depends on the topic and the notation style of the book. For this type of thing, the symbolic stuff shouldn't be OTT.

Try one out in a local library (though most public libraries seem to dreadful for anything above A Level - fortunately I have access to my Uni physics library )

[wossname]

Derby Central Library is the worst in the known universe.

Apart from its reasonable math section it has absolutely nothing else of any consequence at all. The computer section consists of about 400 books on classic ASP, Macromedia Flash, and Photoimpression but nothing about OOP programming, computational logic or ASM.

For saying its in the middle of a major city, its a crap library.