利用Matlab进行极值统计的一个例子——GEV方法

利用Matlab进行极值统计的一个例子——GEV方法科研菜鸟数据和例子均来自于S.Coles,Anintroductiontostatisticalmodelingofextremevalues,Springer,2007一、数据数据是PortPirie,SouthAustralia这个地方的每年海平面的最大值（1923-1987），原始数据可以从R语言包中调出：二、GEV的极大似然估计1、年度极值分布的似然估计广义极值分布为：1/()exp1zGzMatlab中涉及这一方面的函数主要有：parmhat=gevfit(X)returnsmaximumlikelihoodestimatesoftheparametersforthegeneralizedextremevalue(GEV)distributiongiventhedatainX.parmhat(1)istheshapeparameter,K,parmhat(2)isthescaleparameter,sigma,andparmhat(3)isthelocationparameter,mu.[parmhat,parmci]=gevfit(X)returns95%confidenceintervalsfortheparameterestimates.1930195019703.63.84.04.24.44.6YearSeaLevel-----------------------------------------------------------程序----------------------------------------------------clear;clc;load('portpirie.txt');[parmhat,parmci]=gevfit(portpirie);bin=min(portpirie):0.05:max(portpirie);fori=1:length(bin)g(i)=sum(portpirie=bin(i))./length(portpirie);endx=parmhat(1).*((bin-parmhat(3))./parmhat(2));y=-parmhat(1).*log(-log(g));plot(x,y,'ko')holdon;plot(-0.25:0.01:0.1,log(1+(-0.25:0.01:0.1)));xlabel('\xi((z-\mu)/\sigma)');ylabel('-\xilog(-log(G(z)))');-------------------------------------------------------结果---------------------------------------------------------该结果与Coles的结果类似，但尺度参数的CI明显比Coles的结果大。还可利用evim包进行极值统计，evim包的下载地址：~rgencay/evim.html，evim包涉及到这一方面的函数为：out=gev(data,blocksz);wheredataistheinputvectorandblockszisthesizeofthesequentialblocksintowhichthedatawillbesplit.Theoutputoutisastructurethathasthefollowingsevenfields:•out.parests:1×3vectorofestimatedparameters:thefirstelementisestimatedξ,thesecondisestimatedσ,andthethirdisestimatedμ•out.funval:valueofthenegativelog-likelihoodfunction•out.terminated:terminationcondition:1ifsuccessfullyterminated,0otherwise•out.details:detailsofthenonlinearminimizationprocessofthenegativelikelihood•out.varcov:estimatedvariance-covariancematrixoftheparameters•out.parses:estimatedstandarddeviationsoftheestimatedparameters•out.data:vectoroftheextremesoftheblocks-----------------------------------------程序----------------------------------clear;clc;load('portpirie.txt');out=gev(portpirie,1);%out.par_ests(1):shapepara.%out.par_ests(2):scalepara.%out.par_ests(3):positionpara.%95%confidenceintervals:ci_plus=out.par_ests+1.96.*out.par_ses';ci_neg=out.par_ests-1.96.*out.par_ses';%out.par_ses(1):stdofpara.-------------------------------------------------------结果---------------------------------------------------------可见evim包和cole估计的尺度参数的CI均比Matlab自带程序算的小。2、回归水平的估计当0时，{1[log(1)]}pzp当0时，logppzy其中，是位置参数，是尺度参数，是形状参数，且{}pPMaxzp.有时定义log(1)pyp，画pz关于1logpy的图。显然，在这种图中，可以较为方便的看出形状参数的符号。下图中的直线就是根据方框中的公式计算的。也可以直接根据样本来画returnlevel图。将数据从大到小排列，即：(1)(2)()Nzzz，()iz所对应的回归间隔为：111TiN，()izT的图如下图圆点所示。下图中的虚线是95%信度区间，计算方法为：1.96(),1.96()ppppzVarzzVarz其中：()TpppVarzzVz211[(1)log(1)1]Tpppppzyyyy，注意到此时方差-协方差矩阵的对角线顺序为：()()()VarVVarVar-----------------------------------------程序----------------------------------functionreturn_level(par,vcov,data)%returenlevelplot%parofGevfittedtodata%par(1):shapepara.%par(2):scalepara.%par(3):positionpara.%vcov:variance-covariancematrix%data:seriesofmaximaN=length(data);%empiricalreturnlevelemp_rl=sort(data);emp_rp=1./(1-((1:N)./(N+1)));emp_yp=-log((1:N)./(N+1));%modelreturnlevelx=0.001:0.001:0.999;mod_yp=-log(x);mod_rl=par(3)-(par(2)./par(1)).*(1-(mod_yp).^(-par(1)));mod_rp=1./(1-x);%varianceofreturnlevel.fori=1:length(x)yp=-log(x(i));dzp=[par(2).*par(1).^(-2).*(1-yp.^(-par(1)))-par(2).*par(1).^(-1).*yp.^(-par(1)).*log(yp);-par(1).^(-1).*(1-yp.^(-par(1)));1];varzp(i)=dzp'*vcov*dzp;clearypdzp;end%95%confidenceintervalposi_ci=mod_rl+1.96.*sqrt(varzp);neg_ci=mod_rl-1.96.*sqrt(varzp);figure(1)semilogx(1./emp_yp,emp_rl,'k.');holdon;plot(1./mod_yp,mod_rl);plot(1./mod_yp,posi_ci,'k--');plot(1./mod_yp,neg_ci,'k--');xlabel('1/y_p');ylabel('returnlevel');figure(2)semilogx(emp_rp,emp_rl,'k.');holdon;plot(mod_rp,mod_rl);plot(mod_rp,posi_ci,'k--');plot(mod_rp,neg_ci,'k--');xlabel('returnperiod');ylabel('returnlevel');3、模型的诊断PP图经验CDF为：()()/(1)iGziN根据拟合的GEV计算的CDF为：1/()ˆ()exp1izGz画()()ˆ()()iiGzGz的图即为PP图。--------------------------------------------------------程序---------------------------------------------------functionprob_plot(par,data)%p-pplot%parofGevfittedtodata%par(1):shapepara.%par(2):scalepara.%par(3):positionpara.%data:seriesofmaximaN=length(data);%estimateofprobabilityp1=(1:N)./(N+1);p2=gevcdf(sort(data),par(1),par(2),par(3));plot(p1,p2,'k.');holdon;plot(0:1,0:1)xlabel('Empirical');ylabel('Model')QQ图样本quantile即为()iz，其对应的累积概率为：/(1)iN，代入GEV分布，可得模型的quantile为：ˆ1log1pizN画()ˆipzz的图即为QQ图。--------------------------------------------------------程