d9r10.bas

VB常用数值算法集1 VB常用数值算法集1

Attribute VB_Name = "Module1"
Public NDATAT, XT(1000), YT(1000), ARR(1000), AA, ABDEVT
Sub MEDFIT(X(), Y(), NDATA, A, B, ABDEV)
    SX = 0#
    SY = 0#
    SXY = 0#
    SXX = 0#
    For J = 1 To NDATA
        XT(J) = X(J)
        YT(J) = Y(J)
        SX = SX + X(J)
        SY = SY + Y(J)
        SXY = SXY + X(J) * Y(J)
        SXX = SXX + X(J) ^ 2
    Next J
    NDATAT = NDATA
    DEL = NDATA * SXX - SX ^ 2
    AA = (SXX * SY - SX * SXY) / DEL
    BB = (NDATA * SXY - SX * SY) / DEL
    CHISQ = 0#
    For J = 1 To NDATA
        CHISQ = CHISQ + (Y(J) - (AA + BB * X(J))) ^ 2
    Next J
    SIGB = Sqr(CHISQ / DEL)
    B1 = BB
    F1 = ROFUNC(B1)
    B2 = BB + Abs(3 * SIGB) * Sgn(F1)
    F2 = ROFUNC(B2)
1   If F1 * F2 > 0# Then
        BB = 2# * B2 - B1
        B1 = B2
        F1 = F2
        B2 = BB
        F2 = ROFUNC(B2)
        GoTo 1
    End If
    SIGB = 0.01 * SIGB
2   If Abs(B2 - B1) > SIGB Then
        BB = 0.5 * (B1 + B2)
        If BB = B1 Or BB = B2 Then GoTo 3
        F = ROFUNC(BB)
        If F * F1 >= 0# Then
            F1 = F
            B1 = BB
        Else
            F2 = F
            B2 = BB
        End If
        GoTo 2
    End If
3   A = AA
    B = BB
    ABDEV = ABDEVT / NDATA
End Sub
Function ROFUNC(B)
    N1 = NDATAT + 1
    NML = N1 / 2
    NMH = N1 - NML
    For J = 1 To NDATAT
        ARR(J) = YT(J) - B * XT(J)
    Next J
    Call SORT(NDATAT, ARR())
    AA = 0.5 * (ARR(NML) + ARR(NMH))
    Sum = 0#
    ABDEVT = 0#
    For J = 1 To NDATAT
        D = YT(J) - (B * XT(J) + AA)
        ABDEVT = ABDEVT + Abs(D)
        Sum = Sum + XT(J) * Sgn(D)
    Next J
    ROFUNC = Sum
End Function
Sub SORT(N, RA())
      L = Int(N / 2) + 1
      IR = N
      Do
          If L > 1 Then
            L = L - 1
            RRA = RA(L)
          Else
            RRA = RA(IR)
            RA(IR) = RA(1)
            IR = IR - 1
            If IR = 1 Then
              RA(1) = RRA
              Exit Sub
            End If
          End If
          I = L
          J = L + L
          While J <= IR
            If J < IR Then
              If RA(J) < RA(J + 1) Then J = J + 1
            End If
            If RRA < RA(J) Then
              RA(I) = RA(J)
              I = J
              J = J + J
            Else
              J = IR + 1
            End If
          Wend
          RA(I) = RRA
      Loop
End Sub