基于MATLAB环境数学模型参数估计

1、2021-10-181第第6 6讲：环境数学模型参数讲：环境数学模型参数估计方法（一）估计方法（一）1、一元、多元线性回归参数估计 (自学) 3、非线性回归参数估计方法 (自学) 2021-10-182一、参数估计重要性一、参数估计重要性l 环境数学模型要实现求解，确定其参数是重要的一步。环境数学模型要实现求解，确定其参数是重要的一步。 l 模型中的参数有一些是通过试验手段获得、或者根据文献上资料获得、模型中的参数有一些是通过试验手段获得、或者根据文献上资料获得、或者根据经验值获得。或者根据经验值获得。l 上述参数在具体的时间、空间条件下，是否真正的适合某一模型，或上述参数在具体的时间、空间条

2、件下，是否真正的适合某一模型，或者说其适合的程度有多高要靠参数估计的方法评价，并且以这些从各者说其适合的程度有多高要靠参数估计的方法评价，并且以这些从各种途径获得的模型参数作为初值，利用参数估计的方法，获得最适合种途径获得的模型参数作为初值，利用参数估计的方法，获得最适合的模型参数是参数估计的最主要任务。的模型参数是参数估计的最主要任务。l 参数估计的方法多是一些数学上的优化方法，如最小二乘优化法、非参数估计的方法多是一些数学上的优化方法，如最小二乘优化法、非线性优化规划法、梯度最优化算法、线性回归、非线性回归、经验公线性优化规划法、梯度最优化算法、线性回归、非线性回归、经验公式等。式等。20

3、21-10-183二、什么是基于回归分析的参数估计？二、什么是基于回归分析的参数估计？l 回归分析是一种数理统计的方法，用以估计变量之间的相关关系，这回归分析是一种数理统计的方法，用以估计变量之间的相关关系，这种相关关系可能是线性的，也可能是非线性的。种相关关系可能是线性的，也可能是非线性的。l 回归分析过程是根据因变量和自变量的大量观测数据，发现其大致规回归分析过程是根据因变量和自变量的大量观测数据，发现其大致规律，然后用一定的线性或非线性模型去拟合这些观测数据，回归分析律，然后用一定的线性或非线性模型去拟合这些观测数据，回归分析得出因变量和自变量之间的确切数量相依关系式。得出因变量和自变量

4、之间的确切数量相依关系式。l 反过来，如果已经知道因变量和自变量之间数学模型的基本结构，如反过来，如果已经知道因变量和自变量之间数学模型的基本结构，如直线型、幂指数型、多项式型，只是不知道模型中的确切参数，则根直线型、幂指数型、多项式型，只是不知道模型中的确切参数，则根据观测值，利用回归分析的方法将模型参数确定出来的过程就是参数据观测值，利用回归分析的方法将模型参数确定出来的过程就是参数估计的过程。估计的过程。l 利用回归分析技术进行参数估计实际上要比完全意义上的回归分析简利用回归分析技术进行参数估计实际上要比完全意义上的回归分析简单一些。单一些。 2021-10-184三、基于一元线性回归的

5、参数估计三、基于一元线性回归的参数估计l数学形式数学形式2021-10-185三、基于一元线性回归的参数估计三、基于一元线性回归的参数估计l数学形式数学形式2021-10-186三、基于一元线性回归的参数估计三、基于一元线性回归的参数估计l应用举例应用举例2021-10-187三、基于一元线性回归的参数估计三、基于一元线性回归的参数估计l求解思路求解思路2021-10-188三、基于一元线性回归的参数估计三、基于一元线性回归的参数估计l简单程序简单程序将上述过程转化为如下的M文件。S=205080100150200250300350;mu=1.18 2.17 2.82 3.03 3.49 3.

6、65 3.97 4.11 4.22;x=1./S;y=1./mu;X=ones(size(x,1),1),x;ab,bint,r,rint,stats=regress(y,X,0.05)mumax=1/ab(1);ks=ab(2)*mumaxl运行结果运行结果ab = 0.200 ；12.916 ；bint= （略去）；r =（略去）；rint =（略去）；stats =0.9995 15412.8412 0.000000000001 %可见相关系数平方接近1.0，p值接近0mumax = 4.97ks =64.272021-10-189四、基于多元线性回归的参数估计四、基于多元线性回归的参数

7、估计l问题描述问题描述2021-10-1810四、基于多元线性回归的参数估计四、基于多元线性回归的参数估计l应用举例应用举例2021-10-1811四、基于多元线性回归的参数估计四、基于多元线性回归的参数估计l求解思路求解思路2021-10-1812四、基于多元线性回归的参数估计四、基于多元线性回归的参数估计l编程实现编程实现M=10000000;u=0.5;A=20;xx=500;%给出已知条件t=1803004806609001140156018002100240030003600;C=141504506246565783933022121476932;y=log(C.*sqrt(t);x

8、1=1./t;x2=t;X=ones(size(t,1),1),x1,x2;%构造因变量自变量矩阵b012,bint,r,rint,stats=regress(y,X,0.05) %多元线性回归T=xx/u;B=b012(3)*(-1),B=(-1)*b012(2)/T2 %观察两种途径求得的B是否相等？A0=exp(b012(1)-2*B*T);disp(由B算Dx，);Dx=u2/(4*B),disp(由A0算Dx，);Dx=(M/(A0*A*sqrt(4*pi)22021-10-1813四、基于多元线性回归的参数估计四、基于多元线性回归的参数估计l运行结果运行结果b012 =1.0e+0

9、03 *0.01239675394754 -1.24846558474483 -0.00000124753011bint =（略去）；r =（略去）；rint =（略去）stats =1.0e+007 *0.00000009999996 1.09793698679399 0B =0.00124753011206B =0.00124846558474通过B计算出Dx，Dx = 50.06145204457062通过A0计算出Dx，Dx = 50.10180080903253stats中的第一个元素（R2）非常接近1.0，说明多元回归效果非常好，而且stats中的第三个元素（p值）远小于0.05，

10、印证了回归效果好的结论 2021-10-1814五、基于多项式回归的参数估计五、基于多项式回归的参数估计l问题描述问题描述2021-10-1815五、基于多项式回归的参数估计五、基于多项式回归的参数估计l应用举例应用举例2021-10-1816五、基于多项式回归的参数估计五、基于多项式回归的参数估计l编程求解编程求解【求解】分别多项式回归和多元回归，看二者是否有差别t=1 2 2 3 5 6 7 8 9 10 11 12; y= 92 140 213 267 128 343 208 131 188 282 203 143; plot(t,y,k-);gtext(原始曲线);hold on; %

11、进行多项式拟合估计，通过图形观察拟合效果进行多项式拟合估计，通过图形观察拟合效果 p4=polyfit(t,y,4);y4=polyval(p4,t);plot(t,y4,-*);gtext(4次多项式);hold on p6=polyfit(t,y,6);y6=polyval(p6,t);plot(t,y6,-o);gtext(6次多项式);hold on p8=polyfit(t,y,8);y8=polyval(p8,t);plot(t,y8,-square);gtext(8次多项式);hold on p9=polyfit(t,y,9);y9=polyval(p9,t);plot(t,y9

12、,-v);gtext(9次多项式); %进行多元回归分析，通过相关系数判断回归效果好坏进行多元回归分析，通过相关系数判断回归效果好坏 pr4,bint4,r4,rint4,stats4=regress(y,t.4,t.3,t.2,t,ones(size(t,1),1)%进行4次多项式回归及效果 pr6,bint6,r6,rint6,stats6=regress(y,t.6,t.5,t.4,t.3,t.2,t,ones(size(t,1),1) pr8,bint8,r8,rint8,stats8=regress(y,t.8,t.7,t.6,t.5,t.4,t.3,t.2,t,ones(size(

13、t,1),1) pr9,bint9,r9,rint9,stats9=regress(y,t.9,t.8,t.7,t.6,t.5,t.4,t.3,t.2,t,ones(size(t,1),1)%比较多项式回归和多元分析结果是否一致，以比较多项式回归和多元分析结果是否一致，以9次多项式拟合和次多项式拟合和9次多项式回归对比次多项式回归对比ppr=p9-pr9 %如果ppr=0，则说明多项式拟合和多元回归方法的结果没有差别2021-10-1817五、基于多项式回归的参数估计五、基于多项式回归的参数估计l运行结果运行结果2021-10-1818六、基于非线性回归的参数估计六、基于非线性回归的参数估计l

14、问题描述问题描述上述讨论的线性回归中的上述讨论的线性回归中的“线性线性”并非指并非指y与与x的关系，而是指的关系，而是指y是系数是系数b0、b1、b2等的线性函数等的线性函数,在实际科研工作中，在实际科研工作中，y与参数之间的非线性关系更为常与参数之间的非线性关系更为常见。见。2021-10-1819六、基于非线性回归的参数估计六、基于非线性回归的参数估计l应用举例应用举例-12021-10-1820六、基于非线性回归的参数估计六、基于非线性回归的参数估计l问题求解问题求解【求解】本例用非线性回归的方法估计参数首先编制M函数文件，描述非线性方程。function mu=c3fun36(k,s)

15、If nargin2; display(输入参数太多);mumax=k(1);ks=k(2);mu=mumax*s./(ks+s);调用非线性回归函数，估计参数。S=20 5080100150200250300 350;mu=1.18 2.17 2.82 3.03 3.49 3.65 3.97 4.11 4.22;k0=1 10;k,R,J=nlinfit(S,mu,c3fun36,k0);%拟合出参数kdisp(前例估算出mumax = 4.97 ks =64.27);disp(本次估计出的参数:);mumax=k(1),ks=k(2)kci=nlparci(k,R,J);%获得参数k的置信

16、区间mupred, muci=nlpredci(c3fun36,S,k,R,J);%利用新建的预测模型和原来的S估算u的预测值plot(S,mu,k-,S,mupred,b-*);legend(-观测值;*预测值); %预测值和观测值对比预测值和观测值对比gtext(底物浓度S);gtext(比增长速率u) 2021-10-1821六、基于非线性回归的参数估计六、基于非线性回归的参数估计l运行结果运行结果2021-10-1822六、基于非线性回归的参数估计六、基于非线性回归的参数估计l应用举例应用举例-22021-10-1823六、基于非线性回归的参数估计六、基于非线性回归的参数估计l问题求解

17、问题求解【求解】上述解析解含有余误差函数，其手工计算一般要通过查表的方法，而MATLAB中提供了余误差函数的求解函数erfc()，可以直接实现其求解。%首先编制描述解析解模型的函数function C=c3fun39(Dx,t)c0=350; %mg/Lx=1000;%mu=0.6;%m/sC=(c0/2)*(erfc(x-u*t)./(2*sqrt(Dx*t)+exp(u*x/Dx)*erfc(x+u*t)./(2*sqrt(Dx*t);%然后调用主要函数，进行参数估算然后调用主要函数，进行参数估算。t=60*3 9 14 2124 2935 3744 5056 60;C=0.000.05

18、6.00 80.01130.95210.31280.20313.59330.27341.11345.43349.00;Dx0=50;%给出参数初始值Dx=nlinfit(t,C,c3fun39,Dx0);disp(估计出的纵向弥散系数);Dx 2021-10-1824六、基于非线性回归的参数估计六、基于非线性回归的参数估计l求解结果求解结果2021-10-1825七、基于非线性最小二乘优化的参数估计七、基于非线性最小二乘优化的参数估计 l问题描述问题描述基于非线性最小二乘优化的算法在基于非线性最小二乘优化的算法在MATLAB中主要有两个函数中主要有两个函数lsqnonlin()和和lsqcur

19、vefit() 。如如x=lsqnonlin(Fun, x0)的用法就是以初始值的用法就是以初始值x0为搜索起点，发现一个最好为搜索起点，发现一个最好的的x值，该值使函数值，该值使函数Fun平方和最小（即达到目标平方和最小（即达到目标minsum Fun(x).2时的时的x）。如果，将目标函数设置为观测值和预测值的差，将目标函数中的参数设）。如果，将目标函数设置为观测值和预测值的差，将目标函数中的参数设置为自变量，同时将置为自变量，同时将x0设置为参数的初始值，那么得到的结果将是使观测值设置为参数的初始值，那么得到的结果将是使观测值和预测值之差的平方和最小情况下的参数，这个参数是最优的参数。和

20、预测值之差的平方和最小情况下的参数，这个参数是最优的参数。（非原（非原函数）函数）k= lsqcurvefit(Fun,k0,xdata,ydata)任务是以初始值任务是以初始值k0为搜索起点，发现为搜索起点，发现一个最好的一个最好的k值，该值使函数值，该值使函数Fun的值和观测值之差的平方和最小（即达到目的值和观测值之差的平方和最小（即达到目标标minsum (Fun(k,xdata)-ydata).2时的时的k），其中），其中xdata表示自变量的表示自变量的观测值，观测值，ydata表示因变量的观测值，表示因变量的观测值，k表示参数变量。表示参数变量。(原函数）原函数）2021-10-1

21、826七、七、基于非线性最小二乘优化的参数估计基于非线性最小二乘优化的参数估计 l应用举例应用举例2021-10-1827七、七、基于非线性最小二乘优化的参数估计基于非线性最小二乘优化的参数估计 l问题求解问题求解【求解】首先编制非线性模型的M函数。function p=c3fun310(k,xdata)D=xdata; %方便理解,利用xdata变量表示自变量管径 k1=k(1);k2=k(2);k3=k(3);p=k1+k2*D.k3;调用lsqcurvefit ()主函数，开始估算。xdata=0.1 0.150.30.55 0.8 11.21.6;%自变量管道直径ydata=30 55

22、82153294452592930;%因变量价格k0=1 1 1;k=lsqcurvefit (c3fun310,k0,xdata,ydata);%拟合出参数kdisp(估计出的参数:);k1=k(1),k2=k(2),k3=k(3)k1 =27.4152, k2 =406.3915,k3 =1.71372021-10-1828七、七、基于非线性最小二乘优化的参数估计基于非线性最小二乘优化的参数估计 l问题描述问题描述2021-10-1829七、七、基于非线性最小二乘优化的参数估计基于非线性最小二乘优化的参数估计 l问题求解问题求解【求解求解】首先编制非线性模型的首先编制非线性模型的M函数。函

23、数。 function Zmin=c3fun311(k)D=0.1 0.15 0.30.55 0.8 11.2 1.6;%自变量管道直径观测值pdata=30 55 82153 294 452 592 930;%因变量价格观测值k1=k(1);k2=k(2);k3=k(3); p=k1+k2*D.k3; %因变量价格预测值Zmin=pdata-p; %目标函数因变量观测值和预测值的差调用调用lsqnonlin ()主函数，开始估算主函数，开始估算。k0=1 1 1;k=lsqnonlin (c3fun311,k0);%拟合出参数kdisp(估计出的参数: ); k1=k(1),k2=k(2),k3=k(3) k1 =27.4152, k2 =406.3915,k3 =1.7137 2021-10-1830X(km)0929