module1.vb

基于vb.net的最小二乘法多项式系数计算类. 很好用,比Excel还精确.

Imports System.Math
Module Module1
    'This module group 6 classes:
    '**************
    'Regression Class (base class for calculate Least Squares Regressions)
    '  properties: 
    '           xVex, yVec: input, output data points x, y
    '           NumData: number of points
    '           aMat: store coefficients of matrix for build "normal equations"
    '           bVec: store constants of matrix for build "normal equations" also store
    '                 coefficients of regression when Regr()is executed 
    '  methods:
    '           Regr() : Overridable method for calculate regression
    '           BuildMat() : build "normal equations"

    '******************
    'PowerRegression Class (inherits from Regression Class)
    ' y= aX^b 
    'Calculate Power regression   
    '   Methods: Regr(): Set each data point to natural logarithm and calculate regression
    '            EvaluateX() for get Y given X
    '            EvaluateY() for get X given Y
    '****************
    'ExpRegression Class (inherits from Regression Class)
    'y = a e^(bx) 
    'Calculate Exponential regression   
    '   Methods: Regr(): Set each data point y(i) to Natural Logarithm and calculate regression
    '   EvaluateX() for get Y given X
    '   EvaluateY() for get X given Y

    '****************
    'PolyRegression Class (inherits from Regression Class)
    'y = Ao +A1x + A2x^2... 
    ' Calculate polynomial regression
    ' Properties: DegreePoly: input dregree of polynomial
    '             ChiSquare: get Chi Square
    '             Getr2: get Coefficient of Determination
    '             StandardError: get Standard Error
    '
    'Methods: Regr()
    '         Lfit(): get polynomial coefficients fit
    '         EvaluateX() for get Y given X
    '         EvaluateY() for get X given Y

    '*******************
    'Gasuss Jordan Class
    'solve equations systems. Used for solve normal equations
    'properties: aMat:  input coefficients of equations and return 
    '                     inverse matrix when execute gaussj() method 
    '            bVec:  input constants of equations
    '                      and return solution when execute gaussj() method 
    '            Numeq:  input number of equations

    'Methods: gaussj(): solve equations systems

    '*******************
    'Bairstow Class
    'Find reals and complex roots of polynomials for Bairstow's Method
    'used for interpolation: get Y value, given X in polynomial regression
    'Properties:
    '     Degree: input degree of polynomial
    '     Poly: input polynomial 
    '     Roots: store roots
    '     NumRoot: get number of roots of polynomial
    '     Tolerance: precision of calculates
    '     Iterations: Number of iterations 
    'Methods: brstow() : find real and complex roots of polynomial


    'Base Class for calculate Least Squares Regressions
    Public Class Regression
        Protected x(100), y(100) As Double
        Protected a(50, 50), b(50) As Double
        Protected n As Integer
        Public Property xVec() As Array
            Get
                Return x
            End Get
            Set(ByVal value As Array)
                x = value
            End Set
        End Property
        Public Property yVec() As Array
            Get
                Return y
            End Get
            Set(ByVal value As Array)
                y = value
            End Set
        End Property
        Public Property NumData() As Integer
            Get
                Return n
            End Get
            Set(ByVal value As Integer)
                n = value
            End Set
        End Property
        Public Property aMat() As Array
            Get
                Return a
            End Get
            Set(ByVal value As Array)
                a = value
            End Set
        End Property
        Public Property bVec() As Array
            Get
                Return b
            End Get
            Set(ByVal value As Array)
                b = value
            End Set
        End Property

        Public Overridable Sub Regr()
            Dim i As Integer
            Dim g As New GaussJordan

            BuildMat()
            g.NumEq = 2
            g.aMat = a
            g.bVec = b
            g.gaussj()
            a = g.aMat()
            b = g.bVec()
        End Sub
        Public Sub BuildMat()
            Dim i, j, k As Integer

            For i = 1 To n
                For j = 1 To n
                    If i + j > 2 Then
                        For k = 1 To n
                            a(i, j) = a(i, j) + x(k) ^ (i + j - 2)
                        Next k
                    Else
                        a(i, j) = n
                    End If
                Next j
                For k = 1 To n
                    If i > 1 Then
                        b(i) = b(i) + y(k) * x(k) ^ (i - 1)
                    Else
                        b(i) = b(i) + y(k)
                    End If
                Next k
            Next i

        End Sub

    End Class
    'Power regression y= aX^b 
    Public Class PowerRegression : Inherits Regression
        Public Overrides Sub Regr()
            Dim i As Integer
            Dim g As New GaussJordan

            For i = 1 To n
                y(i) = Log(y(i))
                x(i) = Log(x(i))
            Next i

            BuildMat()

            g.NumEq = 2
            g.aMat = a
            g.bVec = b
            g.gaussj()
            a = g.aMat()
            b = g.bVec()
        End Sub
        Public Sub EvaluateX()
            Dim x, y As Double
            Dim xStr As String

            xStr = InputBox("Interpolation (enter x value)")
            x = CDbl(xStr)

            y = Exp(b(1)) * x ^ b(2)
            MsgBox("For x= " & Format(x, "0.000") & vbCrLf & "y= " & Format(y, "0.000"))
        End Sub
        Public Sub EvaluateY()
            Dim x, y As Double
            Dim xStr As String

            xStr = InputBox("Interpolation (enter y value)")
            y = CDbl(xStr)

            x = (y / Exp(b(1))) ^ (1 / b(2))
            MsgBox("For y= " & Format(y, "0.000") & vbCrLf & "x= " & Format(x, "0.000"))
        End Sub

    End Class
    'Exponential regression y = a e^(bx) 
    Public Class ExpRegression : Inherits Regression

        Public Overrides Sub Regr()
            Dim i As Integer
            Dim g As New GaussJordan

            For i = 1 To NumData
                y(i) = Log(y(i))
            Next i

            BuildMat()

            g.NumEq = 2
            g.aMat = a
            g.bVec = b
            g.gaussj()
            a = g.aMat()
            b = g.bVec()
        End Sub
        Public Sub EvaluateX()
            Dim x, y As Double
            Dim xStr As String

            xStr = InputBox("Interpolation (enter x value)")
            x = CDbl(xStr)

            y = Exp(b(1)) * Exp(x * b(2))
            MsgBox("For x= " & Format(x, "0.000") & vbCrLf & "y= " & Format(y, "0.000"))
        End Sub
        Public Sub EvaluateY()
            Dim x, y As Double
            Dim xStr As String

            xStr = InputBox("Interpolation (enter y value)")
            y = CDbl(xStr)

            x = Log(y / Exp(b(1))) / b(2)
            MsgBox("For y= " & Format(y, "0.000") & vbCrLf & "x= " & Format(x, "0.000"))
        End Sub
    End Class
    'Polynomial regression Ao +A1x + A2x^2... 
    Public Class PolyRegression : Inherits Regression
        Private Degree As Integer
        Private Chisq As Double
        Private r2 As Double
        Private StdErr As Double

        Public Property DegreePoly() As Double
            Get
                Return Degree
            End Get
            Set(ByVal value As Double)
                Degree = value
            End Set

        End Property
        Public ReadOnly Property ChiSquare() As Boolean
            Get
                Return Chisq
            End Get
        End Property
        Public ReadOnly Property StandardError() As Boolean
            Get
                Return StdErr
            End Get
        End Property
        Public ReadOnly Property Getr2() As Boolean
            Get
                Return r2
            End Get
        End Property
        Public Overrides Sub Regr()
            Dim i, j As Integer
            Dim covar(20, 20) As Double
            Dim FitTerm(20) As Integer
            Dim xsig(20) As Double
            Dim xStr As String
            Const spread As Integer = 1

            For i = 1 To NumData
                xsig(i) = spread
            Next i
            For i = 1 To Degree + 1
                FitTerm(i) = 1
            Next i
            Lfit(xVec, yVec, xsig, NumData, bVec, FitTerm, Degree, covar, Chisq, r2, StdErr)

        End Sub
        Public Sub Lfit(ByVal x() As Double, ByVal y() As Double, _
                   ByVal xsig() As Double, ByVal ndata As Integer, _
                   ByRef b() As Double, ByVal FitTerm() As Integer, _
                   ByVal Degree As Integer, ByRef covar(,) As Double, _
                   ByRef chisq As Double, ByRef r2 As Double, ByRef StdErr As Double)
            'Given a set of data points x(1..ndata), y(1..ndata) with individual standard deviations
            'sig(1..ndata), use chi2 minimization to fit for some or all of the coefficients b(1..Degree+1) of
            'a function that depends linearly on b, y = SUMi bi AFUNCi(X). The input array FitTerm(1..Degree+1)
            'indicates by nonzero entries those components of a that should be fitted for, and by zero entries
            'those components that should be held fixed at their input values. The program returns values
            'for b(1..Degree+1), chi2 (chisq), coeff of determination (r2), Standard Error (StdErr)
            ' and the covariance matrix covar(l..Degree+1, l..Degree+1).
            '(Parameters held fixed will return zero covariances.)  The routine funcs (x, afunc,Degree+1)
            'returns the Degree+1 basis of polynomial evaluated at x in the array afunc(1..Degree+1)
            'http://www.library.cornell.edu/nr/bookcpdf/c15-1.pdf

            Dim k, kk, j, i, h, m As Integer
            Dim Sumy2, z, v, p, q As Double
            Dim ym, wt, suma, sig2i As Double
            Dim NumTerms, mfit As Integer
            Dim bx(20), afunc(20), bxOrig(20), covarOrig(20, 20) As Double
            Dim g As New GaussJordan
            NumTerms = Degree + 1

            mfit = 0
            For j = 1 To NumTerms
                If FitTerm(j) Then mfit = mfit + 1
            Next j
            If mfit = 0 Then
                MsgBox("lfit: no parameters to be fitted")
                Exit Sub
            End If
            For j = 1 To NumTerms
                For k = 1 To NumTerms
                    covar(j, k) = 0.0
                Next k
                bx(j) = 0.0
            Next j
            For i = 1 To ndata
                funcs(x(i), afunc, NumTerms)
                ym = y(i)
                If mfit < NumTerms Then
                    For j = 1 To NumTerms
                        If Not FitTerm(j) Then ym = ym - b(j) * afunc(j)
                    Next j
                End If
                sig2i = 1.0 / xsig(i) ^ 2

                j = 0
                For h = 1 To NumTerms
                    If FitTerm(h) Then
                        wt = afunc(h) * sig2i
                        j = j + 1
                        k = 0
                        For m = 1 To h
                            k = k + 1
                            If FitTerm(m) Then covar(j, k) = covar(j, k) + wt * afunc(m)
                        Next m
                        bx(j) = bx(j) + ym * wt
                    End If
                Next h
            Next i

            For j = 2 To mfit
                For k = 1 To j - 1
                    covar(k, j) = covar(j, k)
                Next k
            Next j
            bxOrig = bx
            covarOrig = covar
            g.NumEq = mfit
            g.aMat = covar
            g.bVec = bx
            g.gaussj()
            covar = g.aMat()
            bx = g.bVec()
            j = 0
            For h = 1 To NumTerms
                If FitTerm(h) Then
                    j = j + 1
                    b(h) = bx(j)
                End If
            Next h
            chisq = 0.0
            Sumy2 = 0
            For i = 1 To ndata
                funcs(x(i), afunc, NumTerms)
                suma = 0.0
                For j = 1 To NumTerms
                    suma = suma + b(j) * afunc(j)
                Next j
                chisq = chisq + ((y(i) - suma) / xsig(i)) ^ 2
                Sumy2 = Sumy2 + y(i) ^ 2
            Next i