尚德功
(河南省人民胜利渠管理局,河南 新乡 453003)
在实际工作中,经常需要根据一些试验数据求出反映其变化规律的曲线,这些曲线可以是二维、三维等等。这里仅就二维曲线的计算机回归进行探讨。这里的二维是指一个量随另一个量的变化而变化的规律;这里所说的曲线是广义上的曲线,包括直线和折线等。
如图1所示小麦生长期与耗水强度的关系,这条曲线可根据观测试验数据在坐标纸上描点并根据其变化趋势画出平滑的曲线,还可以用插值法、及其它方法来粗略地求出其函数关系式。但如何利用计算机来精确地模拟出曲线,进而准确地求出其函数关系式呢?
图1
对于某些变化规律呈线性或近似线性的函数关系,可用最小二乘法来作线性回归,得出的是一次直线函数关系。对于某些参量之间呈现的不是线性关系,而接近初等几何所包括的几种曲线,如:抛物线、双曲线、对数曲线,指数曲线、椭圆,圆弧等时,也可根据点状图呈现的规律,将曲线假设为相应初等几何曲线的标准方程,来求待定系数。
但是,在科研工作中,常见两维参量之间的对应关系不能用线性和上述初等几何曲线来描述。对于这些函数关系又如何用函数表达式来描述呢?能否找出一种统一的方法,无论是线性还是某种初等几何曲线,或是非初等几何曲线,只要是二维参量函数关系都能根据其试验数据求出其对应的函数表达式呢?进而又能在计算机上得到统一解决呢?
数学中有这样一个定理:若Y=f(X)在点X0的某一邻域内连续且有直到(n+1)阶的连续导数,对X0某个邻域内的任一点X都有:
f(X)=a0+a1(X-X0)+a2(X-X0)2+…+an(X-X0)n+… X∈(X0-δ,X0+δ)
当X0=0时,上式化为:
对任一 X∈[0,δ),显然 X∈(-δ,+δ)。 因此,该定理在[0,δ)上成立。
对于每一组试验数据(X,Y),上式都是关于系数 a0,a1,a2,…,an,…的方程。
即:Y=a0+a1X+a2X2+…+anXn+…
这里,a0,a1,a2,…,an,…皆为一次,问题又转化为求线性方程待定系数问题,即可用最小二乘法来求解,并可用F—检验来观察其拟合情况,察看其显著性水平。(有关用最小二乘法求解方程待定系数的问题本文不再赘述)
将上述二维曲线回归的方法编程在计算机上运行,可以使问题轻松得到解决,可以用C语言来编写,先将二维曲线在X=0点附近的展开式假设为:
f(X)=a0+a1X+a2X2+…+anXn
试验数据代入后,上式成为关于 a0,a1,a2,…,an的待定系数方程,用矩阵最小二乘法求解,得出回归曲线,并显示其相关系数R,(这里0≤R≤l,R值愈接近于0说明回归效果越好;愈接近于1回归效果越差)。可根据R值与0和1的接近程度来评价曲线的回归教果,决定是否增加展开式项数n的值,作下一次回归试验。若回归效果较好,即可进一步实施F一检验。在检验过程中,计算机先后要求输入在0.1,0.05,0.01置信度下的F一检验值。并与E=(N-K)U/KQ的值相比较,判断在相应显著性水平下,曲线的回归效果是否显著。
如果当函数展开式项数n的值增加很大,以至于计算机内存发生溢出,仍不能求出拟合效果良好的曲线时,可根据情况确定在非零点求出该曲线的展开式。这时只要根据程序运行要求,结合实际情况,确定输入X0的值,即可同上运行,直到求出理想的曲线多项式为止。
在科技攻关项目:翟坡试区作物水环境自动监测及预测预报研究课题中,购置的土壤墒情传感器的频率与土壤墒情变化曲线,是生产厂家在实验室里作的,与安装现场实际情况差距很大。为保证所测墒情的准确性,我们对购置的墒情传感器逐个进行了曲线律定。这些曲线就是用该软件进行回归的。
在各传感器的埋设现场取土样烘干得到的土壤含水率数据,与相应的传感器频率观测数据资料如表1。
在应用该软件求其中某一个墒情传感器的变化曲线时,首先将相应的实测数据资料按顺序输入程序的数值语句中,根据实际情况选定适当的X0点(一般情况下先选定X0=0),在X0点附近将函数展开成多顶式。然后,确定展开式的顶数,一般从一项开始试算,每次增加一项。试算期间可以只观测相关系数R的值与0和1的接近程度,不作F-检验。通过试算可以找出R值最接近于0的一次试算,确定其展开式的顶数n的值。再根据n的值进一步运行软件求出曲线,并作F一检验,若能通过F-检验,则该曲线就是所需求的土壤墒情与传感器频率变化函数。若不能通过F-检验,则需要更换展开点X0的值继续试算。
用此方法和该软件回归曲线时,还应注意的一个问题是:数据资料的选取应稍超出实际应用曲线定义域的范围,防止回归出来的曲线在定义域边界附近出现变异现象,以保证所应用曲线的正确性和完整性。
根据上述实际观测资料描出数据点状图,依据点状图发展趋势,补充边界数据如表2。
表2 补充边界数据表
将上述资料输入该软件的数据语句,运行后分别得出回归结果如下:
09#(33# 分站20cm埋深)
Y=-0.11+12.22(X-0)-1.48(X-0)2+7.18E-02(X-0)3-1.56E-03(X-0)4+1.26E-05(X-0)5
R=Qe/Syy=0.4792
Et=(n-k)U/KQe)=6.0863
当 α=0.1 时,Fa(5,28)=2.06,Et>Fa(5,28),OK!
当 α=0.05 时,Fa(5,28)=2.56,Et>Fa(5,28),OK!
当 α=0.01 时,Fa(5,28)=3.75,Et>Fa(5,28),OK!
11#(33#分站50cm埋深)
Y=-9.3948E-02+13.0948(X-0)-1.8511(X-0)2+0.11(X-0)3-2.9682E-03(X-0)4+2.9681E-05(X-0)5
R=Qe/Syy=0.4431
Et=(n-k)U/KQe)=5.7825
当 α=0.1 时,Fa(5,23)=2.13,Et>Fa(5,23),OK!
当 α=0.05 时,Fa(5,23)=2.64,Et>Fa(5,23),OK!
当 α=0.01 时,Fa(5,23)=3.94,Et>Fa(5,23),OK!
08#(33#分站80cm埋深)
Y=4.4937E-02+12.1704(X-0)-1.5222(X-0)2+7.9693E-02(X-0)3-1.873E-03(X-0)4+1.6311E-05(X-0)5
R=Qe/Syy=0.4886
Et=(n-k)U/KQe)=4.8149
当 α=0.1 时,Fa(5,23)=2.13,Et>Fa(5,23),OK!
当 α=0.05 时,Fa(5,23)=2.64,Et>Fa(5,23),OK!
当 α=0.01 时,Fa(5,23)=3.94,Et>Fa(5,23),OK!
由上述各自的F-捡验可知:上述三条曲线均在α=0.01显著性水平下回归效果显著,墒情传感器能够反映土壤墒情的变化规律,是要求的函数表达式。
为验证所求出的土壤墒情与传感器频率变化曲线,我们继续进行取土样烘干试验,对烘干得出的土壤含水率值与相应的变化函数求出的土壤含水率值进行了比较。如表3。
33#分站20cm埋深09#传感器频率与土壤含水率变化函数
Y=-0.11+12.22(X-0)-1.48(X-0)2+7.18E-02(X-0)3-1.56E-03(X-0)4+1.26E-05(X-0)5
表3 计算值与墒情烘干值比较表
从表3可以看出:虽然回归曲线函数计算值与土样烘干值绝对误差小于0.5的仅占45.45%,但是差值小于l的却占81.82%,对±壤墒情来说,这是可以接受的。因此,该曲线函数能反映土壤墒情传感器频率值与土壤含水率之间的变化规律。
33#分站50cm理深11#传感器频率与土壤含水率变化函数
Y=-9.3948E-02+13.0948(X-0)-1.8511(X-0)2+0.11(X-0)3-2.9682E-03(X-0)4+2.9681E-05(X-0)5
表4 计算值与墒情烘干值比较表
从表4可以看出:回归曲线函数计算值与土样烘干值绝对误差小于0.5的占63.64%,绝对差值小于l的占81.82%,说明用该软件回归出来的曲线函数能较好地反映土壤墒情传感器频率值与土壤含水率之间的变化规律。
33#分站80cm埋深08#传感器频率与土壤含水率变化函数
Y=4.4937E-02+12.1704(X-0)-1.5222(X-0)2+7.9693E-02(X-0)3-1.873E-03(X-0)4+1.6311E-05(X-0)5
表5 计算值与墒情烘干值比较表
从表5可以看出:回归曲线函数计算值与土样烘干值绝对误差小于1的占72.73%,说明用该软件回归出来的曲线函数能较好地反映土壤墒情传感器频率值与土壤含水率之间的变化规律。
通过上述验证,说明用二维参量曲线回归方法和软件得出的曲线函数能较好地反映传感器频率与土壤墒情的变化规律。
用该软件作曲线回归时,只需将原始数据资料中的个别明显非点删除,即可输入计算机进行回归运算,避免了人为干预的因素,曲线的走向完全根据数据分布规律而定。不象插值法和其它方法那样,先要人为地删除许多所谓不规律的点,再根据观察保留有规律的曲线附近的点,这样难免存在观察者厚此薄彼的现象,以致影响理论曲线函数与客观规律的吻合程度。
[1]马玉芳;陈建华;郝扬满.基于C语言对三次样条函数的求解及程序.价值工程.2011.11
[2]苏俊杰.基于矩阵样条函数的二阶矩阵微分方程数值解法研究.合肥工业大学.2010.03
[3]张晓丹;邵帅;刘钦圣.基于样条函数的光滑支持向量机模型.北京科技大学学报.2012.7
[4]邱慧敏;杨济安.样条函数与样条小波.重庆邮电学院学报(自然科学版).2000-03