form1.frm

计算多元自回归

VERSION 5.00
Object = "{67397AA1-7FB1-11D0-B148-00A0C922E820}#6.0#0"; "MSADODC.OCX"
Object = "{CDE57A40-8B86-11D0-B3C6-00A0C90AEA82}#1.0#0"; "MSDATGRD.OCX"
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   5355
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   8775
   LinkTopic       =   "Form1"
   ScaleHeight     =   5355
   ScaleWidth      =   8775
   StartUpPosition =   3  '窗口缺省
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   495
      Left            =   5880
      TabIndex        =   1
      Top             =   4560
      Width           =   1215
   End
   Begin MSDataGridLib.DataGrid DataGrid1 
      Height          =   2775
      Left            =   600
      TabIndex        =   0
      Top             =   240
      Width           =   2175
      _ExtentX        =   3836
      _ExtentY        =   4895
      _Version        =   393216
      HeadLines       =   1
      RowHeight       =   15
      BeginProperty HeadFont {0BE35203-8F91-11CE-9DE3-00AA004BB851} 
         Name            =   "宋体"
         Size            =   9
         Charset         =   134
         Weight          =   400
         Underline       =   0   'False
         Italic          =   0   'False
         Strikethrough   =   0   'False
      EndProperty
      BeginProperty Font {0BE35203-8F91-11CE-9DE3-00AA004BB851} 
         Name            =   "宋体"
         Size            =   9
         Charset         =   134
         Weight          =   400
         Underline       =   0   'False
         Italic          =   0   'False
         Strikethrough   =   0   'False
      EndProperty
      ColumnCount     =   2
      BeginProperty Column00 
         DataField       =   ""
         Caption         =   ""
         BeginProperty DataFormat {6D835690-900B-11D0-9484-00A0C91110ED} 
            Type            =   0
            Format          =   ""
            HaveTrueFalseNull=   0
            FirstDayOfWeek  =   0
            FirstWeekOfYear =   0
            LCID            =   2052
            SubFormatType   =   0
         EndProperty
      EndProperty
      BeginProperty Column01 
         DataField       =   ""
         Caption         =   ""
         BeginProperty DataFormat {6D835690-900B-11D0-9484-00A0C91110ED} 
            Type            =   0
            Format          =   ""
            HaveTrueFalseNull=   0
            FirstDayOfWeek  =   0
            FirstWeekOfYear =   0
            LCID            =   2052
            SubFormatType   =   0
         EndProperty
      EndProperty
      SplitCount      =   1
      BeginProperty Split0 
         BeginProperty Column00 
         EndProperty
         BeginProperty Column01 
         EndProperty
      EndProperty
   End
   Begin MSAdodcLib.Adodc Adodc1 
      Height          =   735
      Left            =   2640
      Top             =   960
      Width           =   1455
      _ExtentX        =   2566
      _ExtentY        =   1296
      ConnectMode     =   0
      CursorLocation  =   3
      IsolationLevel  =   -1
      ConnectionTimeout=   15
      CommandTimeout  =   30
      CursorType      =   3
      LockType        =   3
      CommandType     =   2
      CursorOptions   =   0
      CacheSize       =   50
      MaxRecords      =   0
      BOFAction       =   0
      EOFAction       =   0
      ConnectStringType=   1
      Appearance      =   1
      BackColor       =   -2147483643
      ForeColor       =   -2147483640
      Orientation     =   0
      Enabled         =   -1
      Connect         =   "Provider=SQLOLEDB.1;Integrated Security=SSPI;Persist Security Info=False;Initial Catalog=wangyan;Data Source=TYUT-60655C3336"
      OLEDBString     =   "Provider=SQLOLEDB.1;Integrated Security=SSPI;Persist Security Info=False;Initial Catalog=wangyan;Data Source=TYUT-60655C3336"
      OLEDBFile       =   ""
      DataSourceName  =   ""
      OtherAttributes =   ""
      UserName        =   ""
      Password        =   ""
      RecordSource    =   "还原量的随机项"
      Caption         =   "Adodc1"
      BeginProperty Font {0BE35203-8F91-11CE-9DE3-00AA004BB851} 
         Name            =   "宋体"
         Size            =   9
         Charset         =   134
         Weight          =   400
         Underline       =   0   'False
         Italic          =   0   'False
         Strikethrough   =   0   'False
      EndProperty
      _Version        =   393216
   End
   Begin VB.Label Label1 
      Height          =   615
      Left            =   120
      TabIndex        =   2
      Top             =   3240
      Width           =   8295
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Private Sub Command1_Click()
Dim zs ' 年径流总量
Dim pj '年径流量平均数
Dim pp
Dim ro1  'ro
Dim ro2
Dim ro3
Dim ro4
Dim fai1 '阶数为1时的回归系数
Dim fai2(1) '阶数为2时的回归系数
Dim fai3(2) '阶数为3时的回归系数
Dim fai4(3)
Dim ybsim(3)   '放ybsim的平方
Dim aic(3) '存放AIC(P)
Dim aic1   '存放最小AIC
Dim h(5) '存放Db(3)低通滤波系数
Dim g(5) '存放Db(3)高通滤波系数
Call dksjk(cn)
  Dim rs As adodb.Recordset
  Set rs = New adodb.Recordset
  rs.Open "select * from 还原量的随机项", cn, adOpenStatic, adLockOptimistic
  Set DataGrid1.DataSource = Adodc1
  n = 0
  While Not rs.EOF
    n = n + 1
  rs.MoveNext
  Wend
 Print n
  ReDim jll(n - 1)
  
  Call dksjk(cn)
  Dim rrs As adodb.Recordset
  Set rrs = New adodb.Recordset
  rrs.Open "select * from 还原量的随机项", cn, adOpenStatic, adLockOptimistic
  Set DataGrid1.DataSource = Adodc1
  i = 0
  zs = 0
  While Not rrs.EOF
  jll(i) = rrs.Fields(0).Value
 ' Print jll(i)
  zs = zs + jll(i)
  i = i + 1
  rrs.MoveNext
  Wend
  pj = zs / n '年径流量平均值
 ' Print pj ' 经测试
For i = 0 To n - 1
pp = pp + (jll(i) - pj) ^ 2
Next i
dat = Sqr(pp / (n - 1)) '年径流量方差
'Print dat  经测试
For i = 1 To n - 1
   pk = pk + (jll(i) - pj) * (jll(i - 1) - pj)
Next i
'Print pk 经测试
ro1 = pk / pp
'Print ro1
fai1 = ro1
For i = 2 To n - 1
   pk1 = pk1 + (jll(i) - pj) * (jll(i - 2) - pj)
Next i

ro2 = pk1 / pp
'Print ro2
fai2(0) = ro1 + ro1 * ro2
fai2(1) = ro1 ^ 2 + ro2
'Print fai2(0)
'Print fai2(1)
For i = 3 To n - 1
 pk2 = pk2 + (jll(i) - pj) * (jll(i - 3) - pj)
Next i
ro3 = pk2 / pp
fai3(2) = (ro3 - ro1 * ro2 + ((ro2 - ro1 ^ 2) * (ro1 - ro1 * ro2) / (ro1 ^ 2 - 1))) / (1 - ro2 ^ 2 + ((ro1 - ro1 * ro2) ^ 2) / (ro1 ^ 2 - 1))
fai3(1) = (ro2 - ro1 ^ 2 - (ro1 - ro1 * ro2) * fai3(2)) / (1 - ro1 ^ 2)
fai3(0) = ro1 - ro2 * fai3(2) - ro1 * fai3(1) '经检验
For i = 4 To n - 1
 pk3 = pk3 + (jll(i) - pj) * (jll(i - 4) - pj)
Next i
ro4 = pk3 / pp
a = 1 - ro2 ^ 2 + (ro1 - ro1 * ro2) ^ 2 / (ro1 ^ 2 - 1)
b = ro1 - ro2 * ro3 + (ro2 - ro3 * ro1) * (ro1 - ro1 * ro2) / (ro1 ^ 2 - 1)
c = ro3 - ro1 * ro2 + (ro2 - ro1 ^ 2) * (ro1 - ro1 * ro2) / (ro1 ^ 2 - 1)
d = ro1 - ro2 * ro3 + (ro2 - ro1 * ro3) * (ro1 - ro1 * ro2) / (ro1 ^ 2 - 1)
e = 1 - ro3 ^ 2 + (ro2 - ro1 * ro3) * (ro2 - ro3 * ro1) / (ro1 ^ 2 - 1)
f = ro4 - ro1 * ro3 + (ro2 - ro1 * ro3) * (ro2 - ro1 ^ 2) / (ro1 ^ 2 - 1)
fai4(3) = (f - d * c / a) / (e - d * b / a)
fai4(2) = c / a - fai4(3) * (b / a)
fai4(1) = (ro2 - ro1 ^ 2 - (ro2 - ro3 * ro1) * fai4(3) - (ro1 - ro1 * ro2) * fai4(2)) / (1 - ro1 ^ 2)
fai4(0) = ro1 - ro3 * fai4(3) - ro2 * fai4(2) - ro1 * fai4(1) '经测试
ybsim(0) = (1 - ro1 * fai1) * dat ^ 2
ybsim(1) = (1 - ro1 * fai2(0) - ro2 * fai2(1)) * dat ^ 2
ybsim(2) = (1 - ro1 * fai3(0) - ro2 * fai3(1) - ro3 * fai3(2)) * dat ^ 2
ybsim(3) = (1 - ro1 * fai4(0) - ro2 * fai4(1) - ro3 * fai4(2) - ro4 * fai4(3)) * dat ^ 2

aic(0) = n * Log(ybsim(0)) + 2 * 1
aic(1) = n * Log(ybsim(1)) + 2 * 2
aic(2) = n * Log(ybsim(2)) + 2 * 3
aic(3) = n * Log(ybsim(3)) + 2 * 4
Debug.Print aic(0), aic(1), aic(2), aic(3)
For i = 0 To 3
Print aic(i)
Next i  '经测试
                        aic1 = aic(3)
                         For i = 0 To 2
                           If aic(i) < aic1 Then
                                 aic1 = aic(i)
                          End If
                         Next i
'Print aic1  经测试
fai1 = Format(fai1, "###.00")
fai2(0) = Format(fai2(0), "###.00")
fai2(1) = Format(fai2(1), "###.00")
fai3(0) = Format(fai3(0), "###.00")
fai3(1) = Format(fai3(1), "###.00")
fai3(2) = Format(fai3(2), "###.00")
fai4(0) = Format(fai4(0), "###.00")
fai4(1) = Format(fai4(1), "###.00")
fai4(2) = Format(fai4(2), "###.00")
fai4(3) = Format(fai4(3), "###.00")
pj = Format(pj, "###.00")
'自回归方程的求解
If aic1 = aic(0) Then
      Label1.Caption = "x(t)=" & pj & "+(" & fai1 & ")*(x(t-1)-" & pj & ")+εt"
ElseIf aic1 = aic(1) Then
     Label1.Caption = "x(t)=" & pj & "+(" & fai2(0) & ")*(x(t-1)-" & pj & ")+(" & fai2(1) & ")*(x(t-2)-" & pj & ")+εt"
ElseIf aic1 = aic(2) Then
     Label1.Caption = "x(t)=" & pj & "+(" & fai3(0) & ")*(x(t-1)-" & pj & ")+(" & fai3(1) & ")*(x(t-2)-" & pj & ")+(" & fai3(2) & ")*(x(t-3)-" & pj & ")+εt"
ElseIf aic1 = aic(3) Then
     Label1.Caption = "x(t)=" & pj & "+(" & fai4(0) & ")*(x(t-1)-" & pj & ")+(" & fai4(1) & ")*(x(t-2)-" & pj & ")+(" & fai4(2) & ")*(x(t-3)-" & pj & ")+(" & fai4(3) & ")*(x(t-4)-" & pj & ")+εt"
End If
Debug.Print ro1, ro2, ro3, ro4
'进行小波分解
h(0) = 0.3327
h(1) = 0.8069
h(2) = 0.4599
h(3) = -0.135
h(4) = -0.0854
h(5) = 0.0352
g(0) = 0.0352
g(1) = 0.0854
g(2) = -0.135
g(3) = -0.4599
g(4) = 0.8069
g(5) = -0.3327
ReDim c0(n - 1)
For i = 0 To n - 1  '为c0,k负值为采样序列
    c0(i) = jll(i)
Next i
l = Int((n - 1) / 2)
ReDim c1(l - 1) '存放c1,k
ReDim d1(l - 1) '存放 d1,k
c1(0) = h(2) * c0(0) + h(3) * c0(1) + h(4) * c0(2) + h(5) * c0(3)
c1(l - 1) = h(0) * c0(2 * (l - 1) - 2) + h(1) * c0(2 * (l - 1) - 1) + h(2) * c0(2 * (l - 1))
c1(l - 2) = h(0) * c0(2 * (l - 2) - 2) + h(1) * c0(2 * (l - 2) - 1) + h(2) * c0(2 * (l - 2)) + h(3) * c0(2 * (l - 2) + 1) + h(4) * c0(2 * (l - 2) + 2)
For i = 1 To l - 2
 c1(i) = c0(2 * i - 2) * h(0) + h(1) * c0(2 * i - 1) + h(2) * c0(2 * i) + h(3) * c0(2 * i + 1) + h(4) * c0(2 * i + 2) + h(5) * c0(2 * i + 3) '边界外补零
Next i
d1(0) = g(2) * c0(0) + g(3) * c0(1) + g(4) * c0(2) + g(5) * c0(3)
d1(l - 1) = g(0) * c0(2 * (l - 1) - 2) + g(1) * c0(2 * (l - 1) - 1) + g(2) * c0(2 * (l - 1))
d1(l - 2) = g(0) * c0(2 * (l - 2) - 2) + g(1) * c0(2 * (l - 2) - 1) + g(2) * c0(2 * (l - 2)) + g(3) * c0(2 * (l - 2) + 1) + g(4) * c0(2 * (l - 2) + 2)
For i = 1 To l - 2
 d1(i) = c0(2 * i - 2) * g(0) + g(1) * c0(2 * i - 1) + g(2) * c0(2 * i) + g(3) * c0(2 * i + 1) + g(4) * c0(2 * i + 2) + g(5) * c0(2 * i + 3) '边界外补零
Next i
'进行小波消噪

End Sub