非线性最小二乘法数据拟合（1）.frm

运用最小二乘法的数值拟合

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   3090
   ClientLeft      =   60
   ClientTop       =   450
   ClientWidth     =   4680
   LinkTopic       =   "Form1"
   ScaleHeight     =   3090
   ScaleWidth      =   4680
   StartUpPosition =   3  '窗口缺省
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   4440
      TabIndex        =   0
      Top             =   7680
      Width           =   1335
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Private Sub Command1_Click()
'PROGRAM D9R8A
    'Driver for routine MRQMIN
    NPT = 100
    MA = 6
    SPREAD = 0.001
    Dim X(100), Y(100), SIG(100), A(6), LISTA(6)
    Dim COVAR(6, 6), ALPHA(6, 6), GUES(6)
    A(1) = 5#: A(2) = 2#: A(3) = 3#: A(4) = 2#: A(5) = 5#: A(6) = 3#
    GUES(1) = 4.5: GUES(2) = 2.2: GUES(3) = 2.8
    GUES(4) = 2.5: GUES(5) = 4.9: GUES(6) = 2.8
    IDUM& = -911
    'First try a sum of two Gaussians
    For I = 1 To 100
        X(I) = 0.1 * I
        Y(I) = 0#
        For J = 1 To 4 Step 3
            Y(I) = Y(I) + A(J) * Exp(-((X(I) - A(J + 1)) / A(J + 2)) ^ 2)
        Next J
        Y(I) = Y(I) * (1# + SPREAD * GASDEV(IDUM&))
        SIG(I) = SPREAD * Y(I)
    Next I
    MFIT = MA
    For I = 1 To MFIT
        LISTA(I) = I
    Next I
    ALAMDA = -1
    For I = 1 To MA
        A(I) = GUES(I)
    Next I
    Call MRQMIN(X(), Y(), SIG(), NPT, A(), MA, LISTA(), MFIT, COVAR(), ALPHA(), MA, CHISQ, ALAMDA)
    K = 1
    ITST = 0
1   Print "Iteration #"; K; "Chi-squared: ";
    Print Format$(CHISQ, "#####.##00");
    Print "   ALAMDA: "; Format$(ALAMDA, ".##E+00")
    Print Tab(5)
    Print Tab(3); "A(1)   A(2)    A(3)    A(4)    A(5)    A(6)"
    For I = 1 To MA
        Print Format$(A(I), "#.###0"),
    Next I
    Print Tab(5)
    K = K + 1
    OCHISQ = CHISQ
    Call MRQMIN(X(), Y(), SIG(), NPT, A(), MA, LISTA(), MFIT, COVAR(), ALPHA(), MA, CHISQ, ALAMDA)
    If CHISQ > OCHISQ Then
        ITST = 0
    ElseIf Abs(OCHISQ - CHISQ) < 0.1 Then
        ITST = ITST + 1
    End If
    If ITST < 2 Then
        GoTo 1
    End If
    ALAMDA = 0#
    Call MRQMIN(X(), Y(), SIG(), NPT, A(), MA, LISTA(), MFIT, COVAR(), ALPHA(), MA, CHISQ, ALAMDA)
    Print Tab(5); "Uncertainties:"
    For I = 1 To MA
        Print Format$(Sqr(COVAR(I, I)), "#.###0"),
    Next I
    Print Tab(5); "Expected results:"
    Print Tab(2); "5.0     2.0     3.0     2.0     5.0     3.0"
End Sub
Static Sub MRQMIN(X(), Y(), SIG(), NDATA, A(), MA, LISTA(), MFIT, COVAR(), ALPHA(), NCA, CHISQ, ALAMDA)
    Dim ATRY(20), BETA(20), DA(20)
    If ALAMDA < 0# Then
        KK = MFIT + 1
        For J = 1 To MA
            IHIT = 0
            For K = 1 To MFIT
                If LISTA(K) = J Then IHIT = IHIT + 1
            Next K
            If IHIT = 0 Then
                LISTA(KK) = J
                KK = KK + 1
            ElseIf IHIT > 1 Then
                Print "Improper permutation in LISTA"
                Exit Sub
            End If
        Next J
        If KK <> MA + 1 Then Print "Improper permutation in LISTA"
        ALAMDA = 0.001
        Call MRQCOF(X(), Y(), SIG(), NDATA, A(), MA, LISTA(), MFIT, ALPHA(), BETA(), NCA, CHISQ)
        OCHISQ = CHISQ
        For J = 1 To MA
            ATRY(J) = A(J)
        Next J
    End If
    For J = 1 To MFIT
        For K = 1 To MFIT
            COVAR(J, K) = ALPHA(J, K)
        Next K
        COVAR(J, J) = ALPHA(J, J) * (1# + ALAMDA)
        DA(J) = BETA(J)
    Next J
    Call GAUSSJ(COVAR(), MFIT, DA())
    If ALAMDA = 0# Then
        Call COVSRT(COVAR(), NCA, MA, LISTA(), MFIT)
        Erase DA, ATRY
        Exit Sub
    End If
    For J = 1 To MFIT
        ATRY(LISTA(J)) = A(LISTA(J)) + DA(J)
    Next J
    Call MRQCOF(X(), Y(), SIG(), NDATA, ATRY(), MA, LISTA(), MFIT, COVAR(), DA(), NCA, CHISQ)
    If CHISQ < OCHISQ Then
        ALAMDA = 0.1 * ALAMDA
        OCHISQ = CHISQ
        For J = 1 To MFIT
            For K = 1 To MFIT
                ALPHA(J, K) = COVAR(J, K)
            Next K
            BETA(J) = DA(J)
            A(LISTA(J)) = ATRY(LISTA(J))
        Next J
    Else
        ALAMDA = 10# * ALAMDA
        CHISQ = OCHISQ
    End If
    Erase DA, ATRY
End Sub
Sub MRQCOF(X(), Y(), SIG(), NDATA, A(), MA, LISTA(), MFIT, ALPHA(), BETA(), NALP, CHISQ)
    Dim DYDA(20)
    For J = 1 To MFIT
        For K = 1 To J
            ALPHA(J, K) = 0#
        Next K
        BETA(J) = 0#
    Next J
    CHISQ = 0#
    For I = 1 To NDATA
        Call FGAUSS(X(I), A(), YMOD, DYDA(), MA)
        SIG2I = 1# / (SIG(I) * SIG(I))
        DY = Y(I) - YMOD
        For J = 1 To MFIT
            WT = DYDA(LISTA(J)) * SIG2I
            For K = 1 To J
                ALPHA(J, K) = ALPHA(J, K) + WT * DYDA(LISTA(K))
            Next K
            BETA(J) = BETA(J) + DY * WT
        Next J
        CHISQ = CHISQ + DY * DY * SIG2I
    Next I
    For J = 2 To MFIT
        For K = 1 To J - 1
            ALPHA(K, J) = ALPHA(J, K)
        Next K
    Next J
End Sub
Sub COVSRT(COVAR(), NCVM, MA, LISTA(), MFIT)
    For J = 1 To MA - 1
        For I = J + 1 To MA
            COVAR(I, J) = 0#
        Next I
    Next J
    For I = 1 To MFIT - 1
        For J = I + 1 To MFIT
            If LISTA(J) > LISTA(I) Then
                COVAR(LISTA(J), LISTA(I)) = COVAR(I, J)
            Else
                COVAR(LISTA(I), LISTA(J)) = COVAR(I, J)
            End If
        Next J
    Next I
    SWAP = COVAR(1, 1)
    For J = 1 To MA
        COVAR(1, J) = COVAR(J, J)
        COVAR(J, J) = 0#
    Next J
    COVAR(LISTA(1), LISTA(1)) = SWAP
    For J = 2 To MFIT
        COVAR(LISTA(J), LISTA(J)) = COVAR(1, J)
    Next J
    For J = 2 To MA
        For I = 1 To J - 1
            COVAR(I, J) = COVAR(J, I)
        Next I
    Next J
End Sub
Sub GAUSSJ(A(), N, B())
    Dim IPIV(50), INDXR(50), INDXC(50)
    For J = 1 To N
        IPIV(J) = 0
    Next J
    For I = 1 To N
        BIG = 0#
        For J = 1 To N
            If IPIV(J) <> 1 Then
                For K = 1 To N
                If IPIV(K) = 0 Then
                    If Abs(A(J, K)) >= BIG Then
                        BIG = Abs(A(J, K))
                        IROW = J
                        ICOL = K
                    End If
                ElseIf IPIV(K) > 1 Then
                    Print "Singular matrix"
                End If
                Next K
            End If
        Next J
        IPIV(ICOL) = IPIV(ICOL) + 1
        If IROW <> ICOL Then
            For L = 1 To N
                DUM = A(IROW, L)
                A(IROW, L) = A(ICOL, L)
                A(ICOL, L) = DUM
            Next L
            DUM = B(IROW)
            B(IROW) = B(ICOL)
            B(ICOL) = DUM
        End If
        INDXR(I) = IROW
        INDXC(I) = ICOL
        If A(ICOL, ICOL) = 0# Then Print "Singular matrix."
        PIVINV = 1# / A(ICOL, ICOL)
        A(ICOL, ICOL) = 1#
        For L = 1 To N
            A(ICOL, L) = A(ICOL, L) * PIVINV
        Next L
        B(ICOL) = B(ICOL) * PIVINV
        For LL = 1 To N
            If LL <> ICOL Then
                DUM = A(LL, ICOL)
                A(LL, ICOL) = 0#
                For L = 1 To N
                    A(LL, L) = A(LL, L) - A(ICOL, L) * DUM
                Next L
                B(LL) = B(LL) - B(ICOL) * DUM
            End If
        Next LL
    Next I
    For L = N To 1 Step -1
        If INDXR(L) <> INDXC(L) Then
            For K = 1 To N
                DUM = A(K, INDXR(L))
                A(K, INDXR(L)) = A(K, INDXC(L))
                A(K, INDXC(L)) = DUM
            Next K
        End If
    Next L
End Sub
Function GASDEV(IDUM&)
    Static ISET, GSET
    If ISET = 0 Then
      Do
        V1 = 2# * RAN1(IDUM&) - 1#
        V2 = 2# * RAN1(IDUM&) - 1#
        R = V1 ^ 2 + V2 ^ 2
      Loop While R >= 1# Or R = 0
      FAC = Sqr(-2# * Log(R) / R)
      GSET = V1 * FAC
      GASDEV = V2 * FAC
      ISET = 1
    Else
      GASDEV = GSET
      ISET = 0
    End If
End Function
Static Function RAN1(IDUM&)
    Dim R(97)
    M1& = 259200: IA1& = 7141: IC1& = 54773: RM1 = 0.0000038580247
    M2& = 134456: IA2& = 8121: IC2& = 28411: RM2 = 0.0000074373773
    M3& = 243000: IA3& = 4561: IC3& = 51349
    If IDUM& < 0 Or IFF = 0 Then
        IFF = 1
        IX1& = (IC1& - IDUM&) Mod M1&
        IX1& = (IA1& * IX1& + IC1&) Mod M1&
        IX2& = IX1& Mod M2&
        IX1& = (IA1& * IX1& + IC1&) Mod M1&
        IX3& = IX1& Mod M3&
        For J = 1 To 97
            IX1& = (IA1& * IX1& + IC1&) Mod M1&
            IX2& = (IA2& * IX2& + IC2&) Mod M2&
            R(J) = (CSng(IX1&) + CSng(IX2&) * RM2) * RM1
        Next J
        IDUM& = 1
    End If
    IX1& = (IA1& * IX1& + IC1&) Mod M1&
    IX2& = (IA2& * IX2& + IC2&) Mod M2&
    IX3& = (IA3& * IX3& + IC3&) Mod M3&
    J = 1 + Int((97 * IX3&) / M3&)
    If J > 97 Or J < 1 Then Print "Abnormal exit": Exit Function
    RAN1 = R(J)
    R(J) = (CSng(IX1&) + CSng(IX2&) * RM2) * RM1
End Function
Sub FGAUSS(X, A(), Y, DYDA(), NA)
    Y = 0#
    For I = 1 To NA - 1 Step 3
        ARG = (X - A(I + 1)) / A(I + 2)
        EX = Exp(-ARG ^ 2)
        FAC = A(I) * EX * 2# * ARG
        Y = Y + A(I) * EX
        DYDA(I) = EX
        DYDA(I + 1) = FAC / A(I + 2)
        DYDA(I + 2) = FAC * ARG / A(I + 2)
    Next I
End Sub