jackie2008 发表于 2009-4-3 15:34

帮忙看下非线性回归分析-Marquardt法的vb程序

Dim i%, j%, k%, newError#, localError#, a#(), ay#(), Y1#(), ypar#(), lError#()
    Dim e#, u#, s#, v#, t%                                          'u为约束探测参数,v为缩放常数,e为精度要求,s=-1,0,1,2
    ReDim a(1 To m, 1 To m) As Double                              'a用来记录线形方程组的系数矩阵
    ReDim ay(1 To m) As Double                                     'ay用来记录线形方程组的右侧矩阵
    ReDim lError(n) As Double                                    'lError用来记录迭代过程中的残差
    ReDim ypar(1 To m, 1 To n) As Double
    ReDim Y1(1 To n) As Double
    ReDim r(1 To m) As Double                                       'r记录每次迭代后的待估参数值
    ReDim h(1 To m) As Double                                       'h为每次迭代过程中的参数值与真值之间的误差
   
    Text3.Text = "y=a(exp(-bt)-exp(-ct))"                            '拟合模型
   
    '**********迭代法求参数**********
    localError = 0
'**********根据所给的待估参数值确定每组观测值下的各个待估参数的偏导数**********
    For i = 1 To n Step 1
      ypar(1, i) = Exp(-r(2) * X(i)) - Exp(-r(3) * X(i))
      ypar(2, i) = -X(i) * r(1) * Exp(-r(2) * X(i))
      ypar(3, i) = -X(i) * r(1) * Exp(-r(3) * X(i))
      Y1(i) = r(1) * (Exp(-r(2) * X(i)) - Exp(-r(3) * X(i)))
      Y1(i) = Y(i) - Y1(i)
      localError = localError + Y1(i) * Y1(i)                        '利用所给待估参数的初始值所求出的初始残差平方和
    Next i
    newError = localError                                                               '使下面能够进入循环
    u = 0
'**********初始化线形方程组的系数矩阵**********
    For i = 1 To m Step 1
    For j = 1 To m Step 1
      a(i, j) = 0
    Next j
    Next i
'**********初始化线形方程组的系数矩阵**********
    Call InitMat(ay)
'**********求线形方程组的相关矩阵**********
    For i = 1 To m Step 1
    For j = 1 To m Step 1
      ay(i) = ay(i) + ypar(i, j) * Y1(j)
      For k = 1 To n Step 1
            a(i, j) = a(i, j) + ypar(i, k) * ypar(j, k)
      Next k
    Next j
    Next i
'**********确定约束探侧参数的初始值,此确定方法仅为其中一种**********
    For i = 1 To m Step 1
      u = u + a(i, i)
    Next i
    u = (u / (n * m)) ^ (1 / 2)
                                                                              
    For i = 1 To m Step 1
      a(i, i) = a(i, i) + u
    Next i
    '**********调用解线形方程组的函数解方程组,从而得到相关参数**********
    Call LEGauss(m, a, ay)
    For i = 1 To m
      h(i) = ay(i)
    Next i
    t = 0
'**********利用之前计算所得的残差平方和和h进入循环**********
    While (Not IsTerm(e, h))                                       '循环条件,控制循环即迭代过程
         s = -1: v = 10
         Do
                t = t + 1                                           '用来记录迭代次数
                u = u * (v ^ s)                                          '迭代过程中约束探侧参数的变化
                For i = 1 To m Step 1
                  r(i) = r(i) + h(i)
                  Picture2.Print r(i);                                                '打印输出迭代过程中的迭代参数值
                Next i
               
'**********在迭代过程中确定每组观测值下的各个待估参数的偏导数**********
                For i = 1 To n Step 1
                  ypar(1, i) = Exp(-r(2) * X(i)) - Exp(-r(3) * X(i))
                  ypar(2, i) = -X(i) * r(1) * Exp(-r(2) * X(i))
                  ypar(3, i) = -X(i) * r(1) * Exp(-r(3) * X(i))
                  Y1(i) = r(1) * (Exp(-r(2) * X(i)) - Exp(-r(3) * X(i)))
                  Y1(i) = Y(i) - Y1(i)
                Next i
               
'**********求迭代过程中的残差平方和**********
                localError = 0
                For i = 1 To n Step 1
                For j = 1 To m Step 1
                  lError(i) = Y1(i) - ypar(j, i) * h(j)
                Next j
                  localError = localError + lError(i) * lError(i)
                Next i
               
'**********确定迭代过程中线形方程组的系数矩阵**********
                For i = 1 To m Step 1
                For j = 1 To m Step 1
                  a(i, j) = 0
                Next j
                Next i
                Call InitMat(ay)
                For i = 1 To m Step 1
                For j = 1 To m Step 1
                  ay(i) = ay(i) + ypar(i, j) * Y1(j)
                  For k = 1 To n Step 1
                        a(i, j) = a(i, j) + ypar(i, k) * ypar(j, k)
                  Next k
                Next j
                Next i
               
                For i = 1 To m Step 1
                  a(i, i) = a(i, i) + u
                Next i
               
'**********调用函数解线形方程组**********
                Call LEGauss(m, a, ay)
               
'**********将解得的线形方程组的解赋给对应参量**********
                For i = 1 To m
                  h(i) = ay(i)
                Next i
                s = s + 1                                 '循环条件
         Loop While localError >= newError                     '循环条件,控制循环即迭代过程
         newError = localError
    Wend

jackie2008 发表于 2009-4-3 15:38

运行时总是在u = u * (v ^ s)的地方出现溢出错误,程序中的 Call InitMat(ay)为数组初始化函数,LEGauss(m, a, ay)为解线形方程组的函数,请各位大虾帮帮忙
页: [1]
查看完整版本: 帮忙看下非线性回归分析-Marquardt法的vb程序