Attribute VB_Name = "modMethod"
'复相关分析
Option Explicit
'求简单相关系数
'x(1 To n)：变量，n为观测次数，已知
'y(1 To n)：变量，n为观测次数，已知
'R：相关系数，计算结果
Public Sub YRelation(x() As Double, y() As Double, R As Double)
Dim Xa As Double, Ya As Double, Sxx As Double, Sxy As Double, Syy 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)
Syy = Syy + (y(I) - Ya) ^ 2
Next I
R = Sxy / Sqr(Sxx * Syy)
End Sub
'全主元高斯-约当消去法求逆矩阵
'A(1 To m, 1 To m)：开始存放欲求逆的矩阵，最终存求逆的结果矩阵，m是自变量个数
Public Sub Invert(a() As Double)
Dim n As Integer, ep As Double
Dim I As Integer, J As Integer, K As Integer
Dim I0 As Integer, J0 As Integer
Dim w As Double, z As Double
Dim b(1 To 100) As Double, c(1 To 100) As Double
Dim p(1 To 100) As Double, Q(1 To 100) As Double
n = UBound(a, 1)
ep = 0.0000000001
For K = 1 To n
w = 0#
For I = K To n
For J = K To n
If Abs(a(I, J)) > Abs(w) Then
w = a(I, J): I0 = I: J0 = J
End If
Next J
Next I
If Abs(w) < ep Then
MsgBox "全主元素的绝对值小于0.0000000001，矩阵是奇异的！"
Exit Sub
End If
If I0 K Then
For J = 1 To n
z = a(I0, J): a(I0, J) = a(K, J): a(K, J) = z
Next J
End If
If J0 K Then
For I = 1 To n
z = a(I, J0): a(I, J0) = a(I, K): a(I, K) = z
Next I
End If
p(K) = I0: Q(K) = J0
For J = 1 To n
If J = K Then
b(J) = 1 / w: c(J) = 1
Else
b(J) = -a(K, J) / w: c(J) = a(J, K)
End If
a(K, J) = 0#: a(J, K) = 0#
Next J
For I = 1 To n
For J = 1 To n
a(I, J) = a(I, J) + c(I) * b(J)
Next J
Next I
Next K
For K = n To 1 Step -1
I0 = p(K): J0 = Q(K)
If I0 K Then
For I = 1 To n
z = a(I, I0): a(I, I0) = a(I, K): a(I, K) = z
Next I
End If
If J0 K Then
For J = 1 To n
z = a(J0, J): a(J0, J) = a(K, J): a(K, J) = z
Next J
End If
Next K
End Sub
'求复相关系数和偏相关系数
'x(1 To n, 1 To m)：自变量，已知。n是观测次数，m是自变量的个数
'xx(1 To n)：当前自变量
'y(1 To n)：因变量，已知
'yy(1 To n)：当前因变量
'a(1 To m, 1 To m)：法方程的系数矩阵，计算结果
'b(0 To m)：回归系数，计算结果
'R：复相关系数，计算结果
'Ry(1 To m)：偏相关系数，计算结果
'Rx(1 To m,1 To m)：自变量之间的简单相关系数，计算结果
'F：F检验值，计算结果
't(1 To m)：偏相关系数的t检验值，计算结果
Public Sub YMulti(x() As Double, xx() As Double, y() As Double, _
yy() As Double, a() As Double, b() As Single, R As Double, _
Ry() As Double, Rx() As Double, F As Double, t() As Double)
Dim I As Integer, J As Integer, K As Integer
Dim Xa(1 To 100) As Double, Ya As Double, S2 As Double
Dim Smm(1 To 100, 1 To 100) As Double, Sy(1 To 100) As Double
Dim RR As Double, VV As Double, WW As Double
Dim Rs(1 To 100) As Double
Dim Syy As Double '离差平方和
Dim Q As Double, U As Double 'Q是残差平方和，U是回归平方和
n = UBound(x, 1) 'N是观测数据的组数，即x的行数
m = UBound(x, 2) 'M是自变量的个数，也是x的列数
For I = 1 To m
For K = 1 To n
Xa(I) = Xa(I) + x(K, I)
Next K
Xa(I) = Xa(I) / n '自变量的平均值
Next I
For K = 1 To n
Ya = Ya + y(K)
Next K
Ya = Ya / n '因变量的平均值
'建立法方程
For I = 1 To m
For J = 1 To m
For K = 1 To n
Smm(I, J) = Smm(I, J) + (x(K, I) - Xa(I)) * (x(K, J) - Xa(J))
Next K
Next J
Next I
For I = 1 To m
For K = 1 To n
Sy(I) = Sy(I) + (x(K, I) - Xa(I)) * (y(K) - Ya)
Next K
Next I
For I = 1 To m
For J = 1 To m
a(I, J) = Smm(I, J) 'a（）是法方程的系数矩阵
Next J
Next I
Invert a '求法方程系数矩阵的逆矩阵
'求多元线性回归系数b（）
For I = 1 To m
For J = 1 To m
b(I) = b(I) + a(I, J) * Sy(J)
Next J
Next I
b(0) = Ya
For I = 1 To m
b(0) = b(0) - b(I) * Xa(I)
Next I
'Syy是离差平方和
For K = 1 To n
Syy = Syy + (y(K) - Ya) ^ 2
Next K
'U是回归平方和
For I = 1 To m
U = U + b(I) * Sy(I)
Next I
'R是复相关系数
R = Sqr(U / Syy)
'Q是残差平方和
Q = Syy - U
S2 = Q / (n - m - 1)
F = (U / m) / S2 'F检验值
'求偏相关系数：
'（1）求出Y与每个自变量的简单相关系数
For I = 1 To m
For J = 1 To n
xx(J) = x(J, I)
Next J
YRelation xx, y, RR
Rs(I) = RR
Next I
'（2）求出自变量之间的简单相关系数
For I = 1 To m
For J = 1 To n
yy(J) = x(J, I)
Next J
For J = 1 To m
For K = 1 To n
xx(K) = x(K, J)
Next K
YRelation xx, yy, RR
Rx(I, J) = RR
Next J
Next I
'（3）求出与Y有关的偏相关系数
For I = 1 To m
VV = 1: WW = 1
For J = 1 To m
If I J Then VV = VV * Rs(J)
If I J Then WW = WW * (1 - Rs(J) ^ 2)
Next J
For J = 1 To m
If I J Then VV = VV * Rx(I, J)
If I J Then WW = WW * (1 - Rx(I, J) ^ 2)
Next J
Ry(I) = (Rs(I) - VV) / Sqr(WW)
Next I
For I = 1 To m
t(I) = Abs(Ry(I) * Sqr(n - 3) / Sqr(Abs(1 - Ry(I) ^ 2)))
Next I
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：上侧概率（ '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
