2环境统计常见数据分析方法的MATLAB实现及应用

环境统计常见数据分析方法的MATLAB实现及应用第二讲一、参数估计方法线性回归非线性回归网格搜索一、参数估计方法——基于线性回归/非线性回归、网格搜索1、线性回归y=X×beta+ε其中y为n×1阶因变量的观测值向量；X为n×p阶自变量组的观测值矩阵；beta为p×1阶系数向量；ε为n×1阶的随机误差向量；n为观测次数；p为影响因素个数。如果p=1则为一元线性回归，若p1则为多元线性回归。MATLAB中调用函数：b=regress(y,X)或[b,bint,r,rint,stats]=regress(y,X,alpha)其中b为估计的系数，bint为b的估计区间；r为回归残差，rint为r的估计区间；向量stats给出依次给出了R2统计量、F值以及P值；上述参数是在置信度为100(1-alpha)情况下得出的（此时p应该小于alpha，模型才成立）。另外，如果回归模型中没有考虑常数项，则上述调用格式中的X为由n×p阶自变量组的观测值构成的矩阵（每一列表示一个因素），如果回归模型中包含常数项，则X为由n×(p+1)阶矩阵，其第一列全部为1，后面p列由自变量组的观测值构成（每列表示一个因素）。一、参数估计方法——基于线性回归/非线性回归、网格搜索线性回归-举例均匀流场中，向河流投放某种示踪剂M=10kg，在下游500m处测得一组浓度Ci和对应时间ti的数据，河流平均流速ux=0.5m/s，河流横断面A=20m2，500m处的示踪剂浓度分布如表。试估计参数Dx。表500m处不同时刻的示踪剂浓度计算结果时间(min)358111519263035405060时间(s)1803004806609001140156018002100240030003600C(mg/m3)141504506246565783933022121476932瞬时投放一维扩散模型方程Ci=ixixixtDtuxtDAM4)(exp42一、参数估计方法——基于线性回归/非线性回归、网格搜索求解思路得到Ciit=ixixxtDtuxDAM4)(exp42令A0=xDAM4，xxDuB42，T=xux，则有：Ciit=iittTBA20)(exp两边取对数，有：ln[Ciit]=iiiiBttBTBTAttTBA/2]ln[)(]ln[2020再次令Yi=ln[Ciit]，b0=lnA0+2BT，b1=－BT2，b2=－B，X1i=1/ti，X2i=ti，则可得Yi=b0+b1X1i+b2X2i对上式根据多元线性回归原理，可估计得到b0、b1、b2，从而估计得到Dx。瞬时投放一维扩散模型方程Ci=ixixixtDtuxtDAM4)(exp42一、参数估计方法——基于线性回归/非线性回归、网格搜索编程实现M=10000000;u=0.5;A=20;xx=500;%给出已知条件t=[1803004806609001140156018002100240030003600]';C=[141504506246565783933022121476932]';y=log(C.*sqrt(t));x1=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是否相等？B=(-1)*b012(2)/T^2%观察两种途径求得的B是否相等？A0=exp(b012(1)-2*B*T);disp('由B算Dx，');Dx=u^2/(4*B),disp('由A0算Dx，');Dx=(M/(A0*A*sqrt(4*pi)))^2一、参数估计方法——基于线性回归/非线性回归、网格搜索2、非线性回归上述讨论的线性回归中的“线性”并非指y与x的关系，而是指y是系数b0、b1、b2等的线性函数,在实际科研工作中，y与参数之间的非线性关系更为常见。非线性模型记作：),...,(),,...,,(,),(2121kmbetabetabetabetaxxxxebetaxfy期中y对回归参数beta是非线性的；e为随机误差，服从正态分布（N(0,σ2)）。MATLAB中实现非线性回归时参数估计时，需要编制描述非线性关系的M函数，其中参数beta一般表示为向量的形式：beta=[beta(1)，beta(2)，…，beta(k)]，自变量x一般也表示为矩阵关系，其中x(:,1)表示第一个自变量，x(:,2)表示第二个自变量，依次等等。一、参数估计方法——基于线性回归/非线性回归、网格搜索非线性回归-MATLAB函数MATLAB的调用函数为：[beta,r,J]=nlinfit(X,y,Fun,beta0)其中beta为估计的参数；r为此参数下的回归残差构成的向量；J为雅可比矩阵；X，y分别为观测的自变量和因变量；Fun为进行回归分析的非线性函数模型，可以用M函数定义，也可以用inline()函数定义；beta0为估计参数的初始值。利用上述输出结果，调用函数nlparci()可以得到非线性模型估计参数的95%置信度下的置信区间，调用格式为：betaci=nlparci(beta,r,J)其中betaci表示非线性模型估计参数beta的95%置信度下的置信区间。同样，利用上述输出结果，调用nlpredci()函数可估计非线性模型响应值（因变量）的置信区间，调用格式：[ypred,yci］=nlpredci(Fun,xinput,beta,r,J)输入参数列表中xinput表示自变量的取值（可以是观测值或者根据需要人为给定值）；输出参数中ypred表示对应xinput的因变量预测值；yci表示ypred在95%置信度下的置信区间半宽度，即其置信区间为[ypred-yci，ypred+yci]。一、参数估计方法——基于线性回归/非线性回归、网格搜索非线性回归-举例-nlinfit已知某种化学物质在环境中的降解速度和多个环境因素或者条件有关，它们之间的定量关系可用非线性模型表示为：34231253211/xbxbxbbxxbyrate。其中b1～b5表示要估计的未知参数；x1～x5表示环境因素；yrate表示降解速度。自变量和因变量观测值如表。序号x1x2x3yrate序号x1x2x3yrate1470300108.55710080652.54228580103.798470190654.3534703001204.82910030054134470801200.02101003001208.5547080102.7511100801200.0561001901014.39122853001011.32实验实测环境因素和反应速度数值序号x1x2x3yrate序号x1x2x3yrate1470300108.55710080652.54228580103.798470190654.3534703001204.82910030054134470801200.02101003001208.5547080102.7511100801200.0561001901014.39122853001011.32一、参数估计方法——基于线性回归/非线性回归、网格搜索【求解】①定义模型的M函数，并给出参数初始值beta0=[b10,b20,b30,b40,b50]，然后调用nlinfit()函数得到估计的参数beta、回归残差r、雅可比矩阵J。②利用以上输出结果以及函数nlparci()得到非线性模型估计参数的95%置信度下的置信区间。③调用nlpredci()函数得到非线性模型响应值（因变量）的置信区间。M函数程序如下：functionyrate=c2fun213(b,x)x1=x(:,1);x2=x(:,2);x3=x(:,3);yrate=(b(1)*x2-x3/b(5))./(1+b(2)*x1+b(3)*x2+b(4)*x3);%数组点运算一、MATLAB基本数学运算X=[47030010285801047030012047080120470801010019010100806547019065100300541003001201008012028530010];%定义自变量xy=[8.55003.79004.82000.02002.750014.39002.54004.350013.00008.50000.050011.3200]';%定义因变量ybeta0=[10.50.20.12];%给出参数初始值[beta,r,J]=nlinfit(X,y,'c2fun213',beta0)%调用函数求取参数betaci=nlparci(beta,r,J)%求参数95%置信度下的估计区间xinput=[47030010;2858010;470300120;47080120;4708010];%给出自变量一些值[ypred,yci]=nlpredci('c2fun213',xinput,beta,r,J)%得到因变量的估计区间运行结果：beta=[1.38710.06890.04550.12201.0874]betaci=-0.75413.5282；-0.03770.1755；-0.03180.1228；-0.06020.3042；-0.61262.7873ypred=[8.43153.99044.95710.01182.6603]yci=[0.24590.22190.16440.16670.1419]二、显著性检验——基于方差分析非线性回归-举例2（自己练习）【例】已知一维连续源模型：)0(,0),()0(,),0()0(,0)0,(022ttCtCtCxxCxCuxCDtCx其解析解为：tDutxerfcDuxtDutxerfcCtxCxx2exp22),(0其中C表示污染物浓度(mg/L)；C0表示初始时刻污染物浓度(mg/L)；u表示流速(m/s)；Dx表示纵向弥散系数（m2/s）；x表示距离（m）；t表示时间(s)。已知某条小河u=0.6m/s，在一次中上游稳定释放某示踪剂，形成稳定初始浓度C0=350.0mg/L，在其下游1000m处监测，得到数据。请据此估算纵向弥散系数Dx。下游1000m处的罗丹明浓度监测值时间/min3914212429353744505660浓度/mg/L0.000.056.0080.01130.95210.31280.20313.59330.27341.11345.43349.00二、显著性检验——基于方差分析非线性回归-举例2（自己练习）【求解】上述解析解含有余误差函数，其手工计算一般要通过查表的方法，而MATLAB中提供了余误差函数的求解函数erfc()，可以直接实现其求解。%首先编制描述解析解模型的函数functionC=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*[3914212429353744505660];C=[0.000.056.0080.01130.95210.31280.20313.59330.27341.11345.43349.00];Dx0=50;%给出参数初始值Dx=nlinfit(t,C,'c3fun39',Dx0);disp('估计出的纵向弥散系数');Dx二、显著性检验——基于方差分析非线性回归-举例2（自己练习）%根据估算出的参数，重新进行