Attribute VB_Name = "modMethod"
'分布假设检验
Option Explicit
'正态分布检验时合并频数小于5的区间
'x0是区间的左端点
'x1是区间的右端点
'y是区间的频数
'L合并的次数
Public Sub Merge(x0() As Double, x1() As Double, _
y() As Double, L As Integer)
Dim n As Integer, I As Integer, J As Integer
n = UBound(x0, 1)
I = 1: L = 0
100:
If I = n - L Then GoTo 200
If y(I) < 5 Then '是否需要合并
y(I) = y(I) + y(I + 1) '合并频数
x1(I) = x1(I + 1) '合并右端点
L = L + 1 '记录一次合并
For J = I + 1 To n - 1 '将区间端点和频数数组向前移位
y(J) = y(J + 1)
x0(J) = x0(J + 1)
x1(J) = x1(J + 1)
Next J
GoTo 100 '转去检查新y(I)是否还小于5
End If
I = I + 1 '对下一个频数进行检查
GoTo 100
200:
If y(n - L) < 5 Then '检查最后一个频数，如果也小于5则向前合并
y(n - L - 1) = y(n - L - 1) + y(n - L)
x1(n - L - 1) = x1(n - L)
L = L + 1
End If
End Sub
'二项式或泊松分布检验时合并频数小于5的区间
'y：区间的观测频数
'y1：区间的理论频数
'L：合并的次数
Public Sub Merge_B(y() As Double, y1() As Double, L As Integer)
Dim n As Integer, I As Integer, J As Integer
n = UBound(x0, 1)
I = 1: L = 0
100:
If I = n - L Then GoTo 200
If y(I) < 5 Then '是否需要合并
y(I) = y(I) + y(I + 1) '合并观测频数
y1(I) = y1(I) + y1(I + 1) '合并理论频数
L = L + 1 '记录一次合并
For J = I + 1 To n - 1 '将区间端点和频数数组向前移位
y(J) = y(J + 1)
y1(J) = y1(J + 1)
Next J
GoTo 100 '转去检查新y(I)是否还小于5
End If
I = I + 1 '对下一个频数进行检查
GoTo 100
200:
If y(n - L) < 5 Then '检查最后一个频数，如果也小于5则向前合并
y(n - L - 1) = y(n - L - 1) + y(n - L)
y1(n - L - 1) = y1(n - L - 1) + y1(n - L)
L = L + 1
End If
End Sub
'计算正态分布函数
'x：正态偏离点
'F：下侧概率
Public Sub Norm(x As Double, F As Double)
Dim y As Double, ER As Double, Q As Double
Dim A
Const a1 As Double = 0.0705230784
Const a2 As Double = 0.0422820123
Const a3 As Double = 0.0092705272
Const a4 As Double = 0.0001520143
Const a5 As Double = 0.0002765672
Const a6 As Double = 0.0000430638
y = 0.707106781187 * Abs(x)
A = a4 + y * (a5 + y * a6)
A = a3 + y * A
A = a2 + y * A
A = a1 + y * A
ER = 1 - (1 + y * A) ^ (-16)
Q = 0.5 * ER
If x < 0 Then F = 0.5 - Q Else F = 0.5 + Q
End Sub
'假定为正态分布时，求理论频数与观测频数的卡方值
'x0(1 To n)：区间段的左端点值
'x1(1 To n)：区间段的右端点值
'x(1 To n)：区间段的左右端点的平均值
'y(1 To n)：区间段的频数
'X2：卡方检验值，计算结果
'L：合并的次数（n-L是合并后的区间数）
Public Sub Norm_X2(x0() As Double, x1() As Double, x() As Double, _
y() As Double, x2 As Double, L As Integer)
Dim Xa As Double, Ya As Double, Xv As Double
Dim n As Integer, I As Integer
Dim z0 As Double, z1 As Double
Dim p As Double, P0 As Double, P1 As Double
Dim Fe As Double, xx As Double, pp, ff
n = UBound(x, 1)
For I = 1 To n
Xa = Xa + x(I) * y(I): Ya = Ya + y(I)
Next I
Xa = Xa / Ya '平均值
For I = 1 To n
Xv = Xv + (x(I) - Xa) ^ 2 * y(I)
Next I
Xv = Sqr(Xv / Ya) '标准差
Merge x0, x1, y, L '对频数小于5的区间进行合并
For I = 1 To n - L
z0 = (x0(I) - Xa) / Xv '标准化左端点
z1 = (x1(I) - Xa) / Xv '标准化右端点
Norm z0, P0 'P0为左端点的下侧概率
Norm z1, P1 'P1为右端点的下侧概率
p = P1 - P0 '区间段的理论频率
pp = pp + p
Fe = p * Ya '区间段的理论频数
ff = ff + Fe
xx = (y(I) - Fe) ^ 2 / Fe '卡方值
x2 = x2 + xx '卡方值的和
Next I
End Sub
'计算符合二项式分布事件的概率
'n：试验次数
'p：一次试验中事件A出现的概率
'x：n次试验中事件A出现的次数
'Bin：事件A出现x次的概率
Public Sub Bino(n As Integer, p As Single, x As Integer, bin As Double)
Dim Q As Single, nn As Double, NI As Double, m As Double
Dim I As Integer, C As Double
Q = 1 - p
If x = 0 Then
bin = Q ^ n
Exit Sub
End If
nn = 0
For I = 1 To n
nn = nn + Log(I)
Next I 'NN = n!
m = 0
For I = 1 To x
m = m + Log(I)
Next I 'x!
NI = 0
For I = 1 To n - x
NI = NI + Log(I)
Next I '(n-x)!
C = nn - m - NI
C = C + Log(p ^ x) + Log(Q ^ (n - x))
bin = Exp(C)
End Sub
'假定为二项式分布时，求理论频数与观测频数的卡方值
'x(1 To n)：事件次数，n为观测组数
'y(1 To n)：观测频数，n为观测组数
'y1(1 To n)：理论频数
'P1：一次试验中，事件出现的概率
'nn：试验次数
'X2：卡方检验值，计算结果
'L：合并的次数（n-L是合并后的区间数）
Public Sub Bino_X2(x() As Double, y() As Double, y1() As Double, _
P1 As Single, nn As Integer, x2 As Double, L As Integer)
Dim I As Integer, J As Integer, n As Integer, z As Double
Dim mm As Integer, xx As Double, zmm As Double
n = UBound(x, 1) 'n是原始区间数
For I = 1 To n
mm = mm + y(I) 'mm是观测频数
Next I
'按合并前数据求出理论频数
For I = 1 To n
J = x(I)
Bino nn, P1, J, z 'z是理论频率
y1(I) = z * mm 'y1是理论频数
Next I
Merge_B y, y1, L '将观测频数和理论频数同时合并
For I = 1 To n - L
xx = (y(I) - y1(I)) ^ 2 / y1(I) '求卡方值
x2 = x2 + xx
Next I
End Sub
'计算符合泊松分布事件的概率
'L：n * p
'x：n次试验中事件A出现的次数
'Poi：事件A出现x次的概率
Public Sub Poisson(Lamda As Double, x As Integer, Poi As Double)
Dim E As Double, x1 As Double, I As Integer
x1 = 0
For I = 1 To x
x1 = x1 + Log(I)
Next I
E = x * Log(Lamda) - Lamda - x1
Poi = Exp(E)
End Sub
'假定为泊松分布时，求理论频数与观测频数的卡方值
'x(1 To n)：事件次数，n为观测组数
'y(1 To n)：观测频数，n为观测组数
'y1(1 To n)：理论频数
'X2：卡方检验值，计算结果
'L：合并的次数（n-L是合并后的区间数）
Public Sub Poi_X2(x() As Double, y() As Double, y1() As Double, _
x2 As Double, L As Integer)
Dim I As Integer, J As Integer, n As Integer, z As Double
Dim mm As Integer, ff As Integer, xx As Double, zmm As Double
Dim Lamda As Double
n = UBound(x, 1)
For I = 1 To n
mm = mm + y(I) 'mm是观测频数
ff = ff + x(I) * y(I) 'x * f
Next I
Lamda = ff / mm
'按合并前数据求出理论频数
For I = 1 To n
J = x(I)
Poisson Lamda, J, z 'z是理论频率
y1(I) = z * mm 'y1是理论频数
Next I
Merge_B y, y1, L '将观测频数和理论频数同时合并
For I = 1 To n - L
xx = (y(I) - y1(I)) ^ 2 / y1(I) '求卡方值
x2 = x2 + xx
Next I
End Sub
'计算X2分布的分位数
'n：自由度
'Q：上侧概率
'xx：分位数
Public Sub PCX2(n As Integer, Q As Double, xx As Double)
Dim I As Integer, x As Double, p As Double, W As Double
Dim x0 As Double, pp As Double, d As Double
If n = 1 Then
PNorm Q / 2, x: xx = x * x
Exit Sub
End If
If n = 2 Then
xx = -2 * Log(Q)
Exit Sub
End If
p = 1 - Q: PNorm Q, x: W = 2 / (9 * n)
x0 = n * (1 - W + x * Sqr(W)) ^ 3
For I = 1 To 30
CX2 n, x0, pp, d
If d = 0 Then
xx = x0
Exit Sub
End If
xx = x0 - (pp - p) / d
If Abs(xx - x0) < 10 - 6 * Abs(xx) Then Exit Sub Else x0 = xx
Next I
End Sub
'求正态分布的分位数
'Q：上侧概率
'x：分位数
Public Sub PNorm(Q As Double, x As Double)
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
'计算卡方分布函数和概率密度
'n：自由度
'x2：卡方值
'F：下侧概率
'd：概率密度
Public Sub CX2(n As Integer, x2 As Double, F As Double, d As Double)
Dim PIS As Double, x As Double, CHS As Double, u As Double
Dim IAI As Integer, pp As Double, N2 As Integer, I As Integer
Const PI As Double = 3.14159265359
If x2 = 0 Then
F = 0: d = 0: Exit Sub
End If
PIS = Sqr(PI)
x = x2 / 2
CHS = Sqr(x2)
If (n \ 2) * 2 = n Then 'n为偶数
u = x * Exp(-x)
F = 1 - Exp(-x)
IAI = 2
Else 'n为奇数
u = Sqr(x) * Exp(-x) / PIS
Norm CHS, pp '调用正态分布函数计算过程
F = 2 * (pp - 0.5)
IAI = 1
End If
If IAI = n Then GoTo LL1 Else N2 = n - 2
For I = IAI To N2 Step 2
F = F - 2 * u / I
u = x2 * u / I
Next I
LL1:
d = u / x2
End Sub
