基于VB的 灰色模型预测 和 线性回归预测

灰色模型预测GM(1,1)与线性回归预测（一元、多元）
新建一个工程，添加一个模块（.bas），两个命令按钮：
窗体代码：
Option Explicit
Private Sub Command1_Click() '灰色模型预测
Dim Data As String
Data = "2.67,3.13,3.25,3.36,3.56,3.72"
GM1_1_Predict Data
End Sub
Private Sub Command2_Click() '线性回归预测
Dim X1 As String, X2 As String, X3 As String, X4 As String, Y As String
X1 = "100.38,99.7,92.3,87.6,87.17,88.3,92.75,100.6,90.05;" '最后要加上分号;
X2 = "53.24,51.5,50.5,52.4,59.6,59.7,65.2,62.4,53.68;" '最后要加上分号;
X3 = "226,250,281,272,194,180,105,115,250;" '最后要加上分号;
Y = "644,640,517,425,385,401,448,599,462" '最后不要加上分号;请注意！！！
Linear_Regression_Predict X1 & X2 & X3 & Y
End Sub
模块代码：
Option Explicit
Private Sub Calculate_1_AGO(X_0() As Double, X_1() As Double) '做一次累加生成1-AGO
Dim i As Long, TempX As Double, K As Long
K = UBound(X_0)
ReDim X_1(K)
For i = 0 To K
TempX = TempX + X_0(i)
X_1(i) = TempX
Next i
End Sub
Private Sub Calculate_Matrix_B(X_1() As Double, B() As Double) '计算数据矩阵B Dim i As Long, K As Long
K = UBound(X_1) - 1
ReDim B(K, 1)
For i = 0 To K
B(i, 0) = -0.5 * (X_1(i) + X_1(i + 1))
B(i, 1) = 1
Next i
End Sub
Private Sub Calculate_Matrix_YN(X_0() As Double, YN() As Double) '计算数据矩阵YN
Dim i As Long, K As Long
K = UBound(X_0) - 1
ReDim YN(K, 0)
For i = 0 To K
YN(i, 0) = X_0(i + 1)
Next i
End Sub
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' 函数名：Matrix_Transpotation
' 功能：计算矩阵的转置transpotation
' 参数：m - Integer型变量，矩阵的行数
' n - Integer型变量，矩阵的列数
' mtxA - Double型m x n二维数组，存放原矩阵
' mtxAT - Double型n x m二维数组，返回转置矩阵
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Private Sub Matrix_Transpotation(mtxA() As Double, mtxAT() As Double) Dim i As Integer, j As Integer
Dim M As Integer, N As Integer
M = UBound(mtxA, 2)
N = UBound(mtxA, 1)
ReDim mtxAT(M, N)
For i = 0 To M
For j = 0 To N
mtxAT(i, j) = mtxA(j, i)
Next j
Next i
End Sub
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' 函数名：Matrix_Multiplication
' 功能：计算矩阵的乘法multiplication
' 参数：m - Integer型变量，相乘的左边矩阵的行数
' n - Integer型变量，相乘的左边矩阵的列数和右边矩阵的行数
' l - Integer型变量，相乘的右边矩阵的列数
' mtxA - Double型m x n二维数组，存放相乘的左边矩阵
' mtxB - Double型n x l二维数组，存放相乘的右边矩阵
' mtxC - Double型m x l二维数组，返回矩阵乘积矩阵
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Private Sub Matrix_Multiplication(mtxA() As Double, mtxB() As Double, mtxC() As Double)
Dim i As Integer, j As Integer, K As Integer
Dim M As Integer, N As Integer, L As Integer
M = UBound(mtxA, 1): N = UBound(mtxB, 1): L = UBound(mtxB, 2)
ReDim mtxC(M, L)
For i = 0 To M
For j = 0 To L
mtxC(i, j) = 0#
For K = 0 To N
mtxC(i, j) = mtxC(i, j) + mtxA(i, K) * mtxB(K, j)
Next K
Next j
Next i
End Sub
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' 函数名：Matrix_Inversion
' 功能：矩阵求逆
' 参数：n - Integer型变量，矩阵的阶数
' mtxA - Double型二维数组，体积为n x n。

存放原矩阵A；返回时存放其逆矩阵A-1。

' 返回值：Boolean型，失败为False，成功为True
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Private Function Matrix_Inversion(mtxA() As Double) As Boolean
' 局部变量
Dim N As Integer
N = UBound(mtxA)
ReDim nIs(N) As Integer, nJs(N) As Integer
Dim i As Integer, j As Integer, K As Integer
Dim d As Double, P As Double
' 全选主元，消元
For K = 0 To N
d = 0#
For i = K To N
For j = K To N
P = Abs(mtxA(i, j))
If (P > d) Then
d = P
nIs(K) = i
nJs(K) = j
End If
Next j
Next i
' 求解失败
If (d + 1# = 1#) Then
Matrix_Inversion = False
Exit Function
End If
If (nIs(K) <> K) Then
For j = 0 To N
P = mtxA(K, j)
mtxA(K, j) = mtxA(nIs(K), j)
mtxA(nIs(K), j) = P
Next j
End If
If (nJs(K) <> K) Then
For i = 0 To N
P = mtxA(i, K)
mtxA(i, K) = mtxA(i, nJs(K))
mtxA(i, nJs(K)) = P
Next i
End If
mtxA(K, K) = 1# / mtxA(K, K)
For j = 0 To N
If (j <> K) Then mtxA(K, j) = mtxA(K, j) * mtxA(K, K)
Next j
For i = 0 To N
If (i <> K) Then
For j = 0 To N
If (j <> K) Then mtxA(i, j) = mtxA(i, j) - mtxA(i, K) * mtxA(K, j)
Next j
End If
Next i
For i = 0 To N
If (i <> K) Then mtxA(i, K) = -mtxA(i, K) * mtxA(K, K)
Next i
Next K
' 调整恢复行列次序
For K = N To 0 Step -1
If (nJs(K) <> K) Then
For j = 0 To N
P = mtxA(K, j)
mtxA(K, j) = mtxA(nJs(K), j)
mtxA(nJs(K), j) = P
Next j
End If
If (nIs(K) <> K) Then
For i = 0 To N
P = mtxA(i, K)
mtxA(i, K) = mtxA(i, nIs(K))
mtxA(i, nIs(K)) = P
Next i
End If
Next K
' 求解成功
Matrix_Inversion = True
End Function
Private Sub Predicted_Value(ByVal X_1_0 As Double, ByVal u_value As Double, ByVal a_value As Double, K As Long, PV() As Double)
Dim i As Long
ReDim PV(K)
For i = 1 To K + 1
PV(i - 1) = (X_1_0 - u_value / a_value) * Exp(-a_value * (i - 1)) * (1 - Exp(a_value))
Next i
PV(0) = X_1_0
End Sub
Private Sub String_to_Array(Data As String, X_0() As Double) 'Data字符串转为X_0 数组,X_0 是原始序列
Dim Predict_Data() As String, K As Long, i As Long
Predict_Data = Split(Data, ",")
K = UBound(Predict_Data)
For i = 0 To K
X_0(i) = Predict_Data(i)
Next i
End Sub
Private Sub Print_Array(Arrays() As Double, Title As String) '打印数组Dim i As Long
Form1.Print vbCrLf & String(25, "-") & Title & String(25, "-") & vbCrLf
For i = 0 To UBound(Arrays)
Form1.Print Format(Arrays(i), "0.0000") & " ";
Next i
Form1.Print
End Sub
Private Sub Print_String(Arrays() As String, Title As String) '打印字符-Relative_Residual_Error-RRE
Dim i As Long
Form1.Print vbCrLf & String(25, "-") & Title & String(25, "-") & vbCrLf
For i = 0 To UBound(Arrays)
Form1.Print Arrays(i) & " ";
Next i
Form1.Print
End Sub
Private Sub Print_One_String(S As Double, Title As String)
Form1.Print vbCrLf & String(25, "-") & Title & String(25, "-") & vbCrLf
Form1.Print Space(25) & Format(S, "0.#####")
End Sub
Private Sub Absolute_Residual_Error(Array1() As Double, Array2() As Double, ARE() As Double) '绝对残差
Dim K As Long, i As Long
K = UBound(Array1)
ReDim ARE(K)
For i = 0 To K
ARE(i) = Format(Abs(Array2(i) - Array1(i)), "0.0000")
Next i
End Sub
Private Sub Relative_Residual_Error(Array1() As Double, ARE() As Double, RRE() As String) '相对残差
Dim K As Long, i As Long
K = UBound(Array1)
For i = 0 To K
RRE(i) = Format(ARE(i) / Array1(i), "0.000%")
Next i
End Sub
Private Sub Relatedness_Test(ARE() As Double, P As Double, R As Double) '计算相关度
Dim i As Long, K As Long, Min As Double, Max As Double, SumR As Double Dim Ri() As Double
K = UBound(ARE)
ReDim Ri(K)
Min = ARE(0): Max = ARE(0)
For i = 0 To K
If ARE(i) < Min Then Min = ARE(i)
If ARE(i) > Max Then Max = ARE(i)
Next i
For i = 0 To K
Ri(i) = (Min + P * Max) / (ARE(i) + P * Max)
SumR = SumR + Ri(i)
Next i
R = Format(SumR / (K + 1), "0.0000")
End Sub
Private Function Array_Mean(Array1() As Double) As Double '计算平均数Dim i As Long, K As Long
K = UBound(Array1)
For i = 0 To K
Array_Mean = Array_Mean + Array1(i)
Next i
Array_Mean = Array_Mean / (K + 1)
End Function
Private Function Mean_Square_Error(Array1() As Double) As Double '计算均方差Dim Average As Double, i As Long, K As Long, Temp As Double
K = UBound(Array1)
Average = Array_Mean(Array1())
For i = 0 To K
Temp = Temp + (Array1(i) - Average) ^ 2
Next i
Mean_Square_Error = Format(Sqr(Temp / (K + 1)), "0.0000")
End Function
Public Sub GM1_1_Predict(Data As String) 'Data是原始序列字符,以逗号","分开Dim X_0() As Double, X_1() As Double, B() As Double, YN() As Double, PV() As Double
Dim BT() As Double, BTB() As Double, BTBBT() As Double, A() As Double, ARE() As Double
Dim RRE() As String, R As Double
String_to_Array Data, X_0 'Data字符串转为X_0 数组,X_0 是原始序列
Print_Array X_0, "原始数据"
Calculate_1_AGO X_0, X_1 '做一次累加生成X_1
Calculate_Matrix_YN X_0, YN '计算矩阵YN
Calculate_Matrix_B X_1, B '计算矩阵B
Matrix_Transpotation B, BT '计算矩阵B的转置BT
Matrix_Multiplication BT, B, BTB '矩阵BT×B=BTB
If Not Matrix_Inversion(BTB) Then '矩阵BTB求逆,求逆后也是BTB MsgBox "求解失败", vbCritical, "提示"
Exit Sub
End If
Matrix_Multiplication BTB, BT, BTBBT '矩阵BTB×BT = BTBBT
Matrix_Multiplication BTBBT, YN, A '矩阵BTBBT×YN = A ,A(1, 0)=u,A(0, 0)=a
Debug.Print "u=" & A(1, 0) & "," & "a=" & A(0, 0)
Predicted_Value X_1(0), A(1, 0), A(0, 0), UBound(X_1), PV '预测
Print_Array PV, "预测数据"
Absolute_Residual_Error X_0, PV, ARE
Print_Array ARE, "绝对残差"
Relative_Residual_Error X_0, ARE, RRE
Print_String RRE, "相对残差"
Relatedness_Test ARE, 0.5, R
Print_One_String R, "关联度"
Print_One_String Format(Array_Mean(X_0), "0.0000"), "原始数据平均数"
Print_One_String Format(Array_Mean(ARE), "0.0000"), "残差平均数"
Print_One_String Mean_Square_Error(X_0), "原始数据均方差S1"
Print_One_String Mean_Square_Error(ARE), "残差均方差S2"
End Sub
'''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''以下是线性回归预测'''''''''''''''''''''''
'''''''''''''''''''''''''''''''''''''''''''''''''
Private Sub LRP_OriginalDate_To_Array(OriginalDate As String, X() As Double, Y() As Double)
Dim ODS() As String, i As Long, j As Long, ODS_N As Long, X_Row As Long, X_Column As Long
ODS = Split(OriginalDate, ";")
ODS_N = UBound(ODS) - 1
ReDim Xn(ODS_N) As String
For i = 0 To ODS_N
Xn(i) = ODS(i) 'X1,X2,,,Xn
Next i
Dim TempY() As String
TempY = Split(ODS(ODS_N + 1), ",")
Dim TY As Long
TY = UBound(TempY)
ReDim Y(TY, 0)
For j = 0 To TY
Y(j, 0) = TempY(j)
Form1.Print Y(j, 0)
Next j
X_Row = TY
X_Column = UBound(ODS)
ReDim X(X_Row, X_Column)
Dim TempX() As String
For i = 0 To (X_Column - 1)
TempX = Split(Xn(i), ",")
For j = 0 To X_Row '9
X(j, i + 1) = TempX(j)
X(j, 0) = 1
Next j
Next i
Test_Print_Array X()
End Sub
Public Sub Linear_Regression_Predict(OriginalDate As String)
Dim X() As Double, Y() As Double, B() As Double, i As Long
Dim XT() As Double, XTX() As Double, XTXXT() As Double
LRP_OriginalDate_To_Array OriginalDate, X, Y
Matrix_Transpotation X, XT '计算矩阵X的转置XT
Matrix_Multiplication XT, X, XTX '矩阵XT×X=XTX
If Not Matrix_Inversion(XTX) Then '矩阵BTB求逆,求逆后也是BTB MsgBox "求解失败", vbCritical, "提示"
Exit Sub
End If
Matrix_Multiplication XTX, XT, XTXXT '矩阵XTX×XT = XTXXT
Matrix_Multiplication XTXXT, Y, B '矩阵XTXXT×Y = B ,A(1, 0)=u,A(0, 0)=a
For i = 0 To UBound(B)
Print_One_String B(i, 0), "回归系数b" & i
Next i
End Sub
Private Sub Test_Print_Array(X() As Double) '显示二维数组
Dim i%, j%
For i = 0 To UBound(X(), 1)
For j = 0 To UBound(X(), 2)
Form1.Print Format(X(i, j), "000.0000") & " ";
Next j
Form1.Print
Next i
End Sub
灰色预测
线性回归预测。