d9r8a.frm

数据拟合源程序

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   7050
   ClientLeft      =   1455
   ClientTop       =   540
   ClientWidth     =   7710
   LinkTopic       =   "Form1"
   ScaleHeight     =   7050
   ScaleWidth      =   7710
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   6240
      TabIndex        =   0
      Top             =   6480
      Width           =   1095
   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 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
   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