直线拟合.frm

由测量点坐标

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   9705
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   11520
   LinkTopic       =   "Form1"
   MinButton       =   0   'False
   ScaleHeight     =   9705
   ScaleWidth      =   11520
   WindowState     =   2  'Maximized
   Begin VB.CommandButton Command6 
      Caption         =   "退出"
      Height          =   495
      Left            =   13920
      TabIndex        =   4
      Top             =   10560
      Width           =   1215
   End
   Begin VB.CommandButton Command5 
      Caption         =   "绝对值偏差最小的直线拟合"
      Height          =   495
      Left            =   3120
      TabIndex        =   3
      Top             =   10560
      Width           =   2415
   End
   Begin VB.CommandButton Command2 
      Caption         =   "清屏"
      Height          =   495
      Left            =   12600
      TabIndex        =   2
      Top             =   10560
      Width           =   1335
   End
   Begin VB.PictureBox Picture1 
      BackColor       =   &H80000009&
      Height          =   10215
      Left            =   120
      ScaleHeight     =   10155
      ScaleWidth      =   14955
      TabIndex        =   1
      Top             =   120
      Width           =   15015
   End
   Begin VB.CommandButton Command1 
      Caption         =   "最小二乘法直线拟合"
      Height          =   495
      Left            =   120
      TabIndex        =   0
      Top             =   10560
      Width           =   2655
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Dim NDATAT As Integer
Dim XT(1000) As Single
Dim YT(1000) As Single
Dim ARR(1000) As Single
Dim AA As Single
Dim ABDEVT As Single

Private Sub Form_Load()
 Picture1.ScaleMode = vbMillimeters
 Picture1.ScaleHeight = 30
 Picture1.ScaleWidth = 25
 Picture1.ScaleTop = -20
 Picture1.ScaleLeft = -5
 Picture1.DrawWidth = 1
 
End Sub

Private Sub Command1_Click()

  Picture1.CurrentX = 0
  Picture1.CurrentY = 0
  Picture1.Line (-20, 0)-(20, 0), QBColor(12)
  Picture1.Line (0, -20)-(0, 20), QBColor(12)


  NPT = 100
  SPREAD = 0.5
  Dim X(100), Y(100), SIG(100)
  IDUM = -1984
  
 
  For I = 1 To NPT
   X(I) = 0.1 * I
   Y(I) = -2# * X(I) + 1# + SPREAD * GASDEV(IDUM)
   SIG(I) = SPREAD
   
  Next I
  
  Picture1.CurrentX = X(0)
  Picture1.CurrentY = Y(0)
  For I = 1 To NPT
   Picture1.Line -(X(I), Y(I)), QBColor(2)
  Next I
  
  MWT = 1
  Call FIT(X(), Y(), NPT, SIG(), MWT, A, B, SIGA, SIGB, CHI2, Q)
  Picture1.Line (-5, A - 5 * B)-(20, A + 20 * B), QBColor(12)
  Picture1.Line (-5, 1 - 5 * (-2))-(20, 1 + 20 * (-2)), QBColor(2)
  
  
End Sub

Private Sub Command2_Click()
  Picture1.Cls
End Sub

Private Sub Command5_Click()
Picture1.CurrentX = 0
  Picture1.CurrentY = 0
  Picture1.Line (-20, 0)-(20, 0), QBColor(12)
  Picture1.Line (0, -20)-(0, 20), QBColor(12)


  NPT = 100
  SPREAD = 1.2
  Dim X(100), Y(100), SIG(100)
  IDUM = -1984

  For I = 1 To NPT
   X(I) = 0.1 * I
   Y(I) = -2# * X(I) + 1# + SPREAD * GASDEV(IDUM)
   SIG(I) = SPREAD
   
  Next I
  
  Picture1.CurrentX = 0
  Picture1.CurrentY = 0
  For I = 1 To NPT
   Picture1.Line -(X(I), Y(I)), QBColor(2)
  Next I
  
  MWT = 1
 
  Call MEDFIT(X(), Y(), NPT, A, B, ABDEV)
  Picture1.Line (-5, A - 5 * B)-(20, A + 20 * B), QBColor(12)
  
  Picture1.Line (-5, 1 - 5 * (-2))-(20, 1 + 20 * (-2)), QBColor(2)
 
End Sub

Private Sub Command6_Click()
End
End Sub

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
     DIGB = 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() As Single)
    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 + 1
              Else
                J = IR + 1
              End If
          Wend
          RA(I) = RRA
    Loop
End Sub
Function GASDEV(IDUM)
 Static IDET, GSET
 If ISET = 0 Then
    Do
      Randomize
      V1 = 2# * Rnd - 1#
      V2 = 2# * Rnd - 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

Sub FIT(X(), Y(), NDATA, SIG(), MWT, A, B, SIGA, SIGB, CHI2, Q)
    SX = 0#
    SY = 0#
    ST2 = 0#
    B = 0#
    If MWT <> 0 Then
       SS = 0#
       For I = 1 To NDATA
           WT = 1# / (SIG(I) ^ 2)
           SS = SS + WT
           SX = SX + X(I) * WT
           SY = SY + Y(I) * WT
       Next I
     Else
       For I = 1 To NDATA
           SX = SX + X(I)
           SY = SY + Y(I)
       Next I
       SS = NDATA
     End If
     SXOSS = SX / SS
     If MWT <> 0 Then
        For I = 1 To NDATA
            T = (X(I) - SXOSS) / SIG(I)
            ST2 = ST2 + T * T
            B = B + T * Y(I) / SIG(I)
        Next I
      Else
        For I = 1 To NDATA
            T = X(I) - SXOSS
            ST2 = ST2 + T * T
            B = B + T * Y(I)
        Next I
      End If
      B = B / ST2
      A = (SY - SX * B) / SS
      SIGA = Sqr((1# + SX * SX / (SS * ST2)) / SS)
      SIGB = Sqr(1# / ST2)
      CHI2 = 0#
      If MWT = 0 Then
         For I = 1 To NDATA
             CHI2 = CHI2 + (Y(I) - A - B * X(I)) ^ 2
         Next I
         Q = 1#
         SIGDAT = Sqr(CHI2 / (NDATA - 2))
         SIGA = SIGA * SIGDAT
         SIGB = SIGB * SIGDAT
       Else
         For I = 1 To NDATA
             CHI2 = CHI2 + ((Y(I) - A - B * X(I)) / SIG(I)) ^ 2
         Next I
         Q = GAMMQ(0.5 * (NDATA - 2), 0.5 * CHI2)
      End If
End Sub

Function GAMMQ(A, X)
    If X < 0# Or A <= 0# Then
      Exit Function
    End If
    If X < A + 1# Then
      Call GSER(GAMSER, A, X, GLN)
      GAMMQ = 1# - GAMSER
    Else
      Call GCF(GAMMCF, A, X, GLN)
      GAMMQ = GAMMCF
    End If
End Function
 
Sub GSER(GAMSER, A, X, GLN)
    ITMAX = 100
    EPS = 0.0000003
    GLN = GAMMLN(A)
    If X <= 0# Then
       If X < 0# Then
         Exit Sub
       End If
       GASER = 0#
       Exit Sub
    End If
    AP = A
    Sum = 1# / A
    DEL = Sum
    For N = 1 To ITMAX
        AP = AP + 1#
        DEL = DEL * X / AP
        Sum = Sum + DEL
        If Abs(DEL) < Abs(Sum) * EPS Then GoTo 111
        Next N
111       GAMSER = Sum * Exp(-X + A * Log(X) - GLN)
End Sub

Sub GCF(GAMMCF, A, X, GLN)
    ITMAX = 100
    EPS = 0.000003
    GLN = GAMMLN(A)
    GOLD = 0#
    A0 = 1#
    A1 = X
    B0 = 0#
    B1 = 1#
    FAC = 1#
    For N = 1 To ITMAX
       AN = N
       ANA = AN - A
       A0 = (A1 + A0 * ANA) * FAC
       B0 = (B1 + B0 * ANA) * FAC
       ANF = AN * FAC
       A1 = X * A0 + ANF * A1
       B1 = X * XB0 + ANF * B1
       If A1 <> 0# Then
          FAC = 1# / A1
          G = B1 * FAC
          If Abs((G - GOLD) / G) < EPS Then GoTo 100
          GOLD = G
        End If
    Next N
100 GAMMCF = Exp(-X + A * Log(X) - GLN) * G
          
End Sub

Function GAMMLN(XX)
    Dim COF(6)
    COF(1) = 76.18009173
    COF(2) = -86.50532033
    COF(3) = 24.01409822
    COF(4) = -1.231739516
    COF(5) = 0.00120858003
    COF(6) = -0.00000536382
    STP = 2.50662827465
    HALF = 0.5
    ONE = 1#
    FPF = 5.5
    X = XX - ONE
    TMP = X - FPF
    TMP = (X + HALF) * Log(TMP) - TMP
    SER = ONE
    For J = 1 To 6
        X = X + ONE
        SER = SER + COF(J) / X
    Next J
    GAMMLN = TMP + Log(STP * SER)
End Function