Polynomials [Resolved]

[yrwyddfa]

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

Only for 2x2 matrices . . . (?)

[azteched]

Nah I think for N*N - but I can have a look tomorrow in the library.

[yrwyddfa]

LU decomposition looks promising although it is only more efficient in very particular circumstances but still LU decomposition can be more efficient than row reduction.

I'll look into it.

[wossname]

Mmmmm sweeet.

[yrwyddfa]

Methinks Wossy and Class-A drugs go hand in hand . . .

[dinosaur_uk]

How do i convert this code to VB.Net?

[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)

You might be better talking to Wossy as he's already done the conversion

[dinosaur_uk]

For Stupid COders like me who find this conversion a mission....

Public Structure POINT
        Friend x As Double
        Friend y As Double
    End Structure
    Public Function Trend(ByVal 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
        Dim equation As String
        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)
        MsgBox(equation)
    End Function
    Public Function MxMultiplyCV(ByVal Matrix1(,) As Double, ByVal 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(ByVal 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(ByVal 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

[zaza]

Hi,

Can I recommend for those who are interested in this sort of thing, amongst other mathematical topics,

http://mathworld.wolfram.com/Matrix.html

and links therein. They're the chaps who make Mathematica, so they know what they're on about when it comes to maths programming.

zaza

[guineu]

Wossname or Yrwyddfa,

I know I'm late for this party, but I still need guidance in the use of the code you posted [much] earlier in this thread. When I peruse the code, I get concerned that my input Data(n).x values might blow up the execution, as I have thousands of points to plot with nominal x values, Date and Time, and wish to experiment with high orders. One of you suggested using a For loop to generate the .x values. Will this work for large numbers of points?

In an attempt to answer my own question, perhaps I should assume that a "Trend" line generator would, by definition, focus on a recent subset of points immediately prior to the Trend point being plotted. (I believe Excel limits me to ~3000 points for Poly trendlines.)

Am I on the right track?

Thanks,
--Guineu

[yrwyddfa]

I have replied to your PM.

[Tariq Hasan]

This is a very helpful topic and forum for me. I would like to know whether this code can be modified to fit any of the type of trend like linear, log, power, polynominal so on.. along with the value of coefficients. I mean, it will automatically identify the type of equation to best fit the data series and give the values of coefficients. I am also interested to get the R^2 value to evaluate how best the curve has fitted over the data. Regards. Tariq

[Shaggy Hiker]

You can figure out the R^2 easily enough, but you should keep in mind that it is a really bad idea to fit arbitrary polynomials to data in search of the best fit. The name for that is data dredging, and it has the particular problem that any result you get will almost certainly be wrong.

When some type of regression analysis is performed, it is necessary for the person to start with a hypothesis, which will take the form of something like, "this pattern has a linear relationship to this variable" or "this pattern has a logarithmic relationship to this variable". You are then testing whether or not that relationship is true, and how much of the variation in the pattern is due to the variable. You can't say, "what relationship does this variable have to the pattern I am seeing?" because there MUST be such a relationship. In fact, you can generate two random sequences and look for the relationship between them. Given a sufficiently complex polynomial, you are guaranteed that there will be some equation that relates the two, even though they were independent random numbers. So, if you search through enough polynomials, you are guaranteed to find a relationship between a variable and a pattern, regardless of whether the two have any causal relationship whatsoever. That's the problem with data dredges: They ALWAYS find something, that something will only have meaning by pure chance.

I wrote a program where you could put in some pattern, such as fish returning to spawn in a river. You could then put in whatever factors you thought might relate to the number of fish returning to spawn. The program would then use a genetic algorithm to figure out the relationship between the variables and the pattern, complete with R^2 (which is really just a measure of the minimum sum of the distances from the data points and the polynomial line). After a fair amount of testing, I could show that anybody who believed in the result of the program was a total fool. Fortunately, the program returned a population of excellent polynomial fits. While the best answer was a trap for fools, the population of good results could be examined to determine which variables were actually worth studying further.