一元线性m2.bas

这是一个有关概率中的回归分析算法,内有多种算法,欢迎大家使用.

Attribute VB_Name = "modMethod"
Option Explicit
'一元线性回归
'x(1 To n)：自变量，n为观测次数，已知
'y(1 To n)：因变量，n为观测次数，已知
'b0：截距，计算结果
'b1：回归系数，计算结果
'F：F检验值，计算结果
't：t检验值，计算结果
Public Sub Linear(x() As Double, y() As Double, b0 As Single, b1 As Single, _
                F As Double, t As Double)
    Dim Xa As Double, Ya As Double, Sxx As Double, Sxy As Double, Syy As Double
    Dim SSR As Double, SSE As Double, Ur As Integer, Ue As Integer
    Dim Syx2 As Double, Sb As Double, Sb2 As Double, Sx As Double
    Dim n As Integer, I As Double
    n = UBound(x, 1)
    For I = 1 To n
        Xa = Xa + x(I): Ya = Ya + y(I)
    Next I
    Xa = Xa / n: Ya = Ya / n                        '平均值
    For I = 1 To n
        Sxx = Sxx + (x(I) - Xa) ^ 2
        Sxy = Sxy + (x(I) - Xa) * (y(I) - Ya)
    Next I
    b1 = Sxy / Sxx                                  '截距
    b0 = Ya - b1 * Xa                               '回归系数
'由回归所导致的方差
    For I = 1 To n
        SSR = SSR + (b0 + b1 * x(I) - Ya) ^ 2
    Next I
    Ur = 1
    SSR = SSR / Ur
'由剩余所导致的方差
    For I = 1 To n
        SSE = SSE + (y(I) - b0 - b1 * x(I)) ^ 2
    Next I
Debug.Print "回归方差：", SSR, "剩余方差：", SSE
'求F值
    Ue = n - 2
Debug.Print "回归方差自由度：", Ur, "剩余方差自由度：", Ue
    SSE = SSE / Ue
    F = SSR / SSE
Debug.Print "平均回归方差：", SSR, "平均剩余方差：", SSE, "F检验值：", F
'求t值
    For I = 1 To n
        Syx2 = Syx2 + (y(I) - (b0 + b1 * x(I))) ^ 2
    Next I
    Syx2 = Syx2 / (n - 2)
    For I = 1 To n
        Sx = Sx + (x(I) - Xa) ^ 2
    Next I
    Sb2 = Syx2 / Sx: Sb = Sqr(Sb2): t = Abs(b1 / Sb)
End Sub

'求正态分布的分位数
'Q：上侧概率
'x：分位数
Public Sub PNorm(Q, x)
    Dim p As Double, y As Double, z As Double
    Dim b0 As Double, b1 As Double, b2 As Double
    Dim b3 As Double, b4 As Double, b5 As Double
    Dim b6 As Double, b7 As Double, b8 As Double
    Dim b9 As Double, b10 As Double, b As Double
    b0 = 1.570796288: b1 = 0.03706987906
    b2 = -0.0008364353589: b3 = -0.0002250947176
    b4 = 0.000006841218299: b5 = 0.000005824238515
    b6 = -0.00000104527497: b7 = 8.360937017E-08
    b8 = -3.231081277E-09: b9 = 3.657763036E-11
    b10 = 6.936233982E-13
    If Q = 0.5 Then
        x = 0: GoTo PN01
    End If
    If Q > 0.5 Then p = 1 - Q Else p = Q
    y = -Log(4 * p * (1 - p))
    b = y * (b9 + y * b10)
    b = y * (b8 + b): b = y * (b7 + b)
    b = y * (b6 + b): b = y * (b5 + b)
    b = y * (b4 + b): b = y * (b3 + b)
    b = y * (b2 + b): b = y * (b1 + b)
    z = y * (b0 + b): x = Sqr(z)
    If Q > 0.5 Then x = -x
PN01:
End Sub

'计算F分布的分布函数
'n1：自由度，已知
'n2：自由度，已知
'F：F值，已知
'p：下侧概率，所求
'd：概率密度，所求
Public Sub F_DIST(n1 As Integer, N2 As Integer, F As Double, _
            p As Double, d As Double)
    Dim x As Double, u As Double, Lu As Double
    Dim IAI As Integer, IBI As Integer, nn1 As Integer, nn2 As Integer
    Dim I As Integer
    Const PI As Double = 3.14159265359
    If F = 0 Then
        p = 0: d = 0: Exit Sub
    End If
    x = n1 * F / (N2 + n1 * F)
    If (n1 \ 2) * 2 = n1 Then
        If (N2 \ 2) * 2 = N2 Then
            u = x * (1 - x): p = x: IAI = 2: IBI = 2
        Else
            u = x * Sqr(1 - x) / 2: p = 1 - Sqr(1 - x): IAI = 2: IBI = 1
        End If
    Else
        If (N2 \ 2) * 2 = N2 Then
            p = Sqr(x): u = p * (1 - x) / 2: IAI = 1: IBI = 2
        Else
            u = Sqr(x * (1 - x)) / PI
            p = 1 - 2 * Atn(Sqr((1 - x) / x)) / PI: IAI = 1: IBI = 1
        End If
    End If
    nn1 = n1 - 2: nn2 = N2 - 2
    If u = 0 Then
        d = u / F
        Exit Sub
    Else
        Lu = Log(u)
    End If
    If IAI = n1 Then GoTo LL1
    For I = IAI To nn1 Step 2
        p = p - 2 * u / I
        Lu = Lu + Log((1 + IBI / I) * x)
        u = Exp(Lu)
    Next I
LL1:
    If IBI = N2 Then
        d = u / F: Exit Sub
    End If
    For I = IBI To nn2 Step 2
        p = p + 2 * u / I
        Lu = Lu + Log((1 + n1 / I) * (1 - x))
        u = Exp(Lu)
    Next I
    d = u / F
End Sub

'计算F分布的分位数
'n1：自由度，已知
'n2：自由度，已知
'Q：上侧概率，已知
'F：分位数，所求
Public Sub PF_DIST(n1 As Integer, N2 As Integer, _
                Q As Double, F As Double)
    Dim DF12 As Double, DF22 As Double, A As Double, b As Double
    Dim A1 As Double, b1 As Double, p As Double, YQ As Double
    Dim E As Double, FO As Double, pp As Double, d As Double
    Dim GA1 As Double, GA2 As Double, GA3 As Double
    Dim K As Integer
    DF12 = n1 / 2: DF22 = N2 / 2
    A = 2 / (9 * n1): A1 = 1 - A
    b = 2 / (9 * N2): b1 = 1 - b
    p = 1 - Q: PNorm Q, YQ
    E = b1 * b1 - b * YQ * YQ
    If E > 0.8 Then
        FO = ((A1 * b1 + YQ * Sqr(A1 * A1 * b + A * E)) / E) ^ 3
    Else
        lnGamma DF12 + DF22, GA1
        lnGamma DF12, GA2
        lnGamma DF22, GA3
        FO = (2 / N2) * (GA1 - GA2 - GA3 + 0.69315 + (DF22 - 1) * Log(N2) _
            - DF22 * Log(n1) - Log(Q))
        FO = Exp(FO)
    End If
    For K = 1 To 30
        F_DIST n1, N2, FO, pp, d
        If d = 0 Then
            F = FO: Exit Sub
        End If
        F = FO - (pp - p) / d
        If Abs(FO - F) < 0.000001 * Abs(F) Then Exit Sub Else FO = F
    Next K
End Sub

'计算GAMMA函数
'x：自变量
'z：GAMMA函数值
Public Sub GAMMA(x As Double, z As Double)
    Dim H As Double, y As Double, y1 As Double
    H = 1: y = x
LL1:
    If y = 2 Then
        z = H
        Exit Sub
    ElseIf y < 2 Then
        H = H / y: y = y + 1: GoTo LL1
    ElseIf y >= 3 Then
        y = y - 1: H = H * y: GoTo LL1
    End If
    y = y - 2
    y1 = y * (0.005159 + y * 0.001606)
    y1 = y * (0.004451 + y1)
    y1 = y * (0.07211 + y1)
    y1 = y * (0.082112 + y1)
    y1 = y * (0.41174 + y1)
    y1 = y * (0.422787 + y1)
    H = H * (0.999999 + y1)
    z = H
End Sub

'求Gamma函数的对数LogGamma(x)
'x：自变量
'G：Gamma函数的对数
Public Sub lnGamma(x As Double, G As Double)
    Dim y As Double, z As Double, A As Double
    Dim b As Double, b1 As Double, n As Integer
    Dim I As Integer
    If x < 8 Then
        y = x + 8: n = -1
    Else
        y = x: n = 1
    End If
    z = 1 / (y * y)
    A = (y - 0.5) * Log(y) - y + 0.9189385
    b1 = (0.0007663452 * z - 0.0005940956) * z
    b1 = (b1 + 0.0007936431) * z
    b1 = (b1 - 0.002777778) * z
    b = (b1 + 0.0833333) / y
    G = A + b
    If n >= 0 Then Exit Sub
    y = y - 1: A = y
    For I = 1 To 7
        A = A * (y - I)
    Next I
    G = G - Log(A)
End Sub

'计算t分布的分布函数
'n：自由度，已知
'T：t值，已知
'pp：下侧概率，所求
'dd：概率密度，所求
Public Sub T_Dist(n As Integer, t As Double, pp As Double, dd As Double)
    Dim Sign As Integer, TT As Double, x As Double
    Dim p As Double, u As Double, GA1 As Double, GA2 As Double
    Dim IBI As Integer, N2 As Integer, I As Integer
    Const PI As Double = 3.14159265359
    If t = 0 Then
        Call GAMMA(n / 2, GA1): Call GAMMA(n / 2 + 0.5, GA2): pp = 0.5
        dd = GA2 / (Sqr(n * PI) * GA1): Exit Sub
    End If
    If t < 0 Then Sign = -1 Else Sign = 1
    TT = t * t: x = TT / (n + TT)
    If (n \ 2) * 2 = n Then                 'n为偶数
        p = Sqr(x): u = p * (1 - x) / 2
        IBI = 2
    Else                                    'n为奇数
        u = Sqr(x * (1 - x)) / PI
        p = 1 - 2 * Atn(Sqr((1 - x) / x)) / PI
        IBI = 1
    End If
    If IBI = n Then GoTo LL1 Else N2 = n - 2
    For I = IBI To N2 Step 2
        p = p + 2 * u / I
        u = u * (1 + I) / I * (1 - x)
    Next I
LL1:
    dd = u / Abs(t)
    pp = 0.5 + Sign * p / 2
End Sub

'求t分布的分位数
'n：自由度，已知
'Q：上侧概率（<=0.5），已知
'T：分位数，所求
Public Sub PT_DIST(n As Integer, Q As Double, t As Double)
    Dim PIS As Double, DFR2 As Double, C As Double
    Dim Q2 As Double, p As Double, YQ As Double, E As Double
    Dim GA1 As Double, GA2 As Double, GA3 As Double
    Dim T0 As Double, pp As Double, d As Double
    Dim K As Integer
    Const PI As Double = 3.14159265359
    PIS = Sqr(PI): DFR2 = n / 2
    If n = 1 Then
        t = Tan(PI * (0.5 - Q)): Exit Sub
    End If
    If n = 2 Then
        If Q > 0.5 Then C = -1 Else C = 1
        Q2 = (1 - 2 * Q) ^ 2
        t = Sqr(2 * Q2 / (1 - Q2)) * C
        Exit Sub
    End If
    p = 1 - Q: PNorm Q, YQ              '正态分布分位数
    E = (1 - 1 / (4 * n)) ^ 2 - YQ * YQ / (2 * n)
    If E > 0.5 Then
        T0 = YQ / Sqr(E)
    Else
        lnGamma DFR2, GA1: lnGamma DFR2 + 0.5, GA2
        GA3 = Exp((GA1 - GA2) / n)
        T0 = Sqr(n) / (PIS * Q * n) ^ (1 / n) / GA3
    End If
    For K = 1 To 30
        T_Dist n, T0, pp, d
        If d = 0 Then
            t = T0: Exit Sub
        End If
        t = T0 - (pp - p) / d
        If Abs(T0 - t) < 0.000001 * Abs(t) Then _
            Exit Sub Else T0 = t
    Next K
End Sub