1. 引言
拟合函数的形式通常都是显式的,或称之为标准拟合模型公式,这种形式下因变量y完全可以用自变量x及待求参数向量b组成的表达式进行描述,通俗的讲就是拟合模型公式等号左边仅有因变量y,而右边则是自变量x和参数向量b组成的表达式且不包含因变量y,如公式-1示;而一般隐函数拟合则指拟合模型公式中的因变量y无法用自变量x和参数向量b去表达,如公式-2示,模型公式等号的两侧都包含因变量y,而且也无法求出因变量y的解析表达式。
标准拟合模型形式:
一般隐函数模型形式:
上述标准拟合形式对1stOpt而言是其最为基本也最为强大的功能之一;而一般隐函数拟合,虽然难度有所增加,但1stOpt仍然可以自动应对,用户也不需要做特别额外的处理工作。但如果模型如下公式-3所示,其中的中间变量c无法通过求出c的解析表达式从而消除c,进而转化为标准的拟合形式,这种情况该如何处理?
其中:
由于公式-4无法求出c的解析式,因此公式-3中的c无法消除,这样就形成了一种特殊的隐函数拟合模式:每一步数据点的计算都含有一个等式约束。这种特殊的拟合问题1stOpt能否处理?如果能又该如何处理?下面以三个实际案例分别展示一般及特殊隐函数拟合的1stOpt解决之道。
2. 案例展示分析2.1 案例-1:一般隐函数问题
该案例是经典的隐式函数,如公式-5,拟合公式等号两侧均含有因变量y,拟合数据见表-1.
表-1案例数据
x |
1,2,3,4,6,8,10 |
y |
2.525,2.788,1.276,0.810,0.450,0.303,0.225 |
![]()
对这种“常规”隐函数拟合,1stOpt的处理方式与一般拟合问题并无任何区别,虽然内部实际求解过程完全不同。求解代码及结果如下。
1stOpt一般隐函数计算代码:
Algorithm= UGO1; MaxIteration= 3500; Parameterp( 4); Functiony=exp(p4*x)/(p1+y+p2*x^p3); Data; x= 1, 2, 3, 4, 6, 8, 10; y= 2. 525, 2. 788, 1. 276, 0. 810, 0. 450, 0. 303, 0. 225;
1stOpt一般隐函数计算结果:
ObjectiveFunction( Min.): 0.117229486390231 SumSquaredError( SSE): 0.117229485382952 RootofMeanSquareError( RMSE): 0.129410468434442 CorrelationCoef. (R): 0.99176794362072 DeterminationCoef. ( DC): 0.982640135961821 ![]()
前述隐函数公式-5实际上是能够求出y的解析式的,比如用Matlab中的“solve”函数可得到y的显示表达式如下:
由此可以直接进行显示函数的拟合计算,相比隐式函数拟合速度会快很多,代码及结果如下。
1stOpt代码
Parameterp( 4); Functiony=sqrt( 4*exp(p4*x)+p1^ 2+p2^ 2*x^( 2*p3)+ 2*p1*p2*x^p3)/ 2-p2*x^p3/ 2-p1/ 2; Data; x= 1, 2, 3, 4, 6, 8, 10; y= 2. 525, 2. 788, 1. 276, 0. 810, 0. 450, 0. 303, 0. 225; 1stOpt计算结果
Sum Squared Error (SSE): 0.117229483047838Root of Mean Square Error (RMSE): 0.129410467145568Correlation Coef. ( R): 0.991768205569815R-Square: 0.98360417357917Adjusted R-Square: 0.975406260368755Determination Coef. (DC): 0.982640136307615Chi -Square: 0.107695483017389F -Statistic: 59.009979697594两种方式的计算结果是非常接近的,证明1stOpt处理一般隐式函数拟合的方法是合理的,结果也是正确的。这一对比验证说明,当隐式函数确实无法转换成显式格式时,1stOpt的计算结果是可信的。
2.2 案例-2:特殊隐函数之一
该案例数据同表-1,拟合公式如下:
其中:x、y分别为自变量和因变量,c为一中间变量,满足如下关系:
如果能将公式-8中的中间变量c用x、p1和p4表示,也即如果能求出c的解析表达式,再代入公式-7,就成为一般的拟合问题并可轻易求解了,但问题是目前还无法由公式-8求出c的解析式,也就无法将公式-7中的c替换掉从而求解,甚至也无法按案例-1的因变量隐式拟合方式去处理,这种情况该如何处理呢?
由公式-8可看出,中间变量c对应于自变量x,因为无法求出c的显式表达式,在此将c值系列设定为待求参数,同时增加等式约束条件:
其中:i=1,2,…n,n为拟合数据点数。
实际运行时将上述n个等式约束合并为一个平方和等式约束,即:
这种带特殊约束的拟合问题快捷模式无法实现,可在编程模式下解决。下面给出Fortran及Pascal编程模式下的求解代码。代码中通过关键字“PassParameter”定义并输出约束误差值f。
1stOpt编程模式Fortran代码
MaxIteration = 2500; Variable x,y;Parameter p(1:4); Parameter c(0:6); PassParameter f;StartProgram [Fortran];Subroutine MainModelinteger ireal(8)tem_Errortem_error=0doi=0, DataLength - 1y(i) = (c(i)+x(i)+p4)/(p1+p2*x(i)**p3)tem_error = tem_error + (x(i)-exp(c(i)*x(i))+(p1*x(i)+c(i))**p4)** 2end dof=tem_errorConstrainedResult=(tem_error = 0) End SubroutineEndProgram;Data;x= 1, 2, 3, 4, 6, 8, 10; y= 2.525, 2.788, 1.276, 0.810, 0.450, 0.303, 0.225; 1stOpt编程模式Pascal代码
MaxIteration =2500; Variable x,y;Parameterp( 1: 4); Parameterc( 0: 6); PassParameter f;StartProgram [Pascal];ProcedureMainModel; var i: integer; tem_Error: double; Begintem_error : =0; fori : =0toDataLength -1do beginy[i] : =(c[i] +x[i] +p4) /(p1 +p2 *x[i] ^p3); tem_error : =tem_error +sqr(x[i] -exp(c[i] *x[i]) +(p1 *x[i] +c[i]) ^p4); end; f : =tem_error; ConstrainedResult : =(tem_error =0); End; EndProgram;Data;x =1, 2, 3, 4, 6, 8, 10; y =2.525, 2.788, 1.276, 0.810, 0.450, 0.303, 0.225; 上述两段代码可得出如下一致结果。
1stOpt计算结果
Sum Squared Error (SSE): 0.0466255192773601Root of Mean Square Error (RMSE): 0.081613653687323Correlation Coef. (R): 0.997366647703129R-Square: 0.994740229950578Adjusted R-Square: 0.992110344925867Determination Coef. (DC): 0.993095485553653F-Statistic: -61.2785103151751![]()
2.3 案例-3:特殊隐函数之二
该案例数据与前述相同,拟合公式如下:
其中:x、y分别为自变量和因变量,c为一中间变量,满足如下关系:
很明显,与案例2相比。除了拟合公式由公式-7变为如下的公式-11,其余均完全一致。最大的不同就是拟合模型公式也是隐式函数,即公式等式两端都含有因变量y,且无法求出y的解析表达式,如果加上如公式-12的隐式约束,那么这种“双重隐式“函数拟合问题该怎么处理?
其实非常简单,参照案例2,采取同样的处理方式:增加对应于因变量y数据量的向量参数yy,约束条件如下:
其中:i=1,2,…n,n为拟合数据点数。
实际运行时将公式-12和13各n个等式约束合并为一个平方和等式约束,即:
具体运行代码及结果如下。
1stOpt编程模式Fortran代码
MaxIteration= 2500; Variablex,y; Parameterp( 1: 4); Parameterc( 0: 6), y_( 0: 6)=[ 0, 3]; PassParameterf; StartProgram[Fortran];Subroutine MainModelinteger ireal(8) tem_Error, tem_ytem_error = 0do i = 0, DataLength - 1tem_y = (c(i)+x(i)+(p4+y_(i))**p5)/(p1+p2*x(i)**p3)y(i) = tem_ytem_error = tem_error + (x(i)-exp(c(i)*x(i))+(p1*x(i)+c(i))**p4)**2+(y_(i)-tem_y)**2end dof = tem_errorConstrainedResult = (tem_error = 0)End SubroutineEndProgram;Data;x=1,2,3,4,6,8,10;y=2.525,2.788,1.276,0.810,0.450,0.303,0.225;1stOpt编程模式Pascal代码
Hardness= 2; // 1, 2,... 5MaxIteration= 2500; Variablex,y; Parameterp( 1: 5); Parameterc( 0: 6), y_( 0: 6)=[ 0, 3]; PassParameterf; StartProgram[Pascal];Procedure MainModel;var i: integer;tem_Error, tem_y: double;Begintem_error := 0;for i := 0 to DataLength - 1 do begintem_y := (c[i]+x[i]+(p4+y_[i])^p5)/(p1+p2*x[i]^p3);y[i] := tem_y;tem_error := tem_error + sqr(x[i]-exp(c[i]*x[i])+(p1*x[i]+c[i])^p4)+sqr(y_[i]-tem_y);end;f := tem_error;ConstrainedResult := (tem_error = 0);End;EndProgram;Data;x=1,2,3,4,6,8,10;y=2.525,2.788,1.276,0.810,0.450,0.303,0.225;上述两段代码可得出如下一致结果。
1stOpt计算结果
Sum Squared Error (SSE): 0.0122532419025211Root of Mean Square Error (RMSE): 0.0418385364100083Correlation Coef. (R): 0.999142102993635R-Square: 0.998284941974543Adjusted R-Square: 0.997427412961814Determination Coef. (DC): 0.998185485394226F-Statistic: -372.43512476707![]()
该案例可以看作案例1和案例2的叠加,虽然难度增加很多,但1stOpt仍然可以应对处理。
3. 小结
通过三个实际案例展示了1stOpt处理一般及特殊隐式函数拟合问题。
1stOpt对特殊隐式函数的处理方式是通过定义额外的参数并增加相应的约束去实现的,优点是可以计算任意类型的复杂隐式函数拟合问题,缺点则是待求未知参数数随拟合数据点数增加而增加,从而导致优化求解难度指数级加大、耗时也更多,因此处理实际隐函数拟合问题时,应优先尝试转换成显示形式。
编辑/范瑞强
审核/范瑞强
复核/范瑞强
本文来源:数学中国
关注返回搜狐,查看更多