主题:[讨论]帮忙看下非线性回归分析-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
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