Printable View
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 . . . .
You could use a genetic algorithm. It's an interesting challenge but very fast for high order polynomials.
[deleted because it was obsolete]
Here's my effort: it appears to work identically to Excel:
Enjoy!Code: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
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".
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.Quote:
Originally posted by yrwyddfa
You can derive the relevant expression from the code posted. Look at the last two loops of the Trend function.
Add this to the bottom of the Trend function.Code: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
Cool.
Could it generate this red curve (or an approximation) from the noisy data?
http://www.vbforums.com/attachment.p...postid=1738801
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?
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.
I would'nt call 4 or 5 high order though :D 20 yes. 30 certainly.Quote:
Originally posted by yrwyddfa
Yup; that looks like quite a high order polynomial.
Quote:
How would I call your Trend() function to process this data?
Lets say Raw(0-999) contains the point data.
VB Code:
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 . . .
Well, it's identical.Quote:
y=ax^5 + bx^4 + cx^3 + dx^2 + ex + f
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,
Is this going to end being sanded?
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 :)
Good stuff :) Let me know how you get on.
OMG!!
Dude this RULES. I can scarcely contain myself :D. I had to take a snap for posterity...
http://www.vbforums.com/attachment.p...postid=1832701
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.
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.
:cool: :cool: :cool: :cool: :cool: :cool: :cool: :cool: :cool:
:eek2: :eek2: :eek2: :eek2: :eek2:
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.
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)
What does ph34r mean?
Sorted out the dot drawing problem and optimised your code for .net...
http://www.vbforums.com/attachment.p...postid=1833608
I don't know anything about matrices (except the kind you sleep on :D) 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 :thumb:
This is really cool :cool:
ph34r is l33t-speak for FEAR!
if that's your EEG, then you are one sick kitty!
:lol: yeah, but it looks sweet though right?
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.Quote:
Originally posted by wossname
Sorted out the dot drawing problem and optimised your code for .net...
http://www.vbforums.com/attachment.p...postid=1833608
I don't know anything about matrices (except the kind you sleep on :D) 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 :thumb:
This is really cool :cool:
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.
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.Quote:
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
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.
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? :confused:Quote:
Originally posted by opus
Cool code, but why do I get a different formula back if i run the function a second time?
I'm using the code as posted (including the Debug.Print for the formula).VB Code:
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
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
Add the Checkmark and the word [RESOLVED] to the subject of the first post in your thread
if your question has been answered satisfactorily.
:lol: I bet you were the milk monitor at school weren't you? :DQuote:
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.
it looked like it needed *something*. next time, I'll just send a PM to have it manually [RESOLVED] with the tick.
.Net Framework graphics suck.
I'm going to need to redesign the drawing routines.
DC, and GDI then . . .
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.
If you want the inverse of a matrix, there's faster ways than Gauss-Jordan - read it up in a computational maths textbook :)
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.
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 :) )
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.