利用Matlab进行极值统计的一个例子GPD方法.docx
- 文档编号:24673312
- 上传时间:2023-05-30
- 格式:DOCX
- 页数:15
- 大小:218.47KB
利用Matlab进行极值统计的一个例子GPD方法.docx
《利用Matlab进行极值统计的一个例子GPD方法.docx》由会员分享,可在线阅读,更多相关《利用Matlab进行极值统计的一个例子GPD方法.docx(15页珍藏版)》请在冰豆网上搜索。
利用Matlab进行极值统计的一个例子GPD方法
利用Matlab进行极值统计的一个例子——GPD方法
科研菜鸟
数据和例子均來口PS・Coles,Anintroductiontostaasticalmodelingofextremevalues,Spnnger,2007
一、数据
数据是south-westEngland这个地方的口降雨量(1914-1962),原始数据町以从R语肓包中调出:
>library(isir.皂v)
>data(rain)
>write.table(rainFf±le=nF:
/study/razn.txtn)
s
o寸
Index
二.gpd的极大似然估计
1、超越量的似然估计
广义Pareto分布为:
(z
P{X-z/ 其中.y>0,l+ 分布和广义极值分布的参数关系为: Nlatlab中涉及这一力面的因数主要何: parmhat二gpfit(X)returnsmaximumlikelihoodestimatesoftheparametersforthen^o-parametergeneralizedPareto(GP)distributiongiventhedatainX.pamihat (1)isrhetadindex(shape)parameter,Kandpamihat (2)isthescaleparameter,sigma・gpfitdoesnotfitathreshold(location)parameter・ [parmhat,parmci]=gpfit(X) returns95%confidencemren^lsfortheparameterestimates. 程序 clear;clc; load('rain.txt'); srain=rain(rain>30); [parmhatparmci]=gpfit(srain-30); bin=min(srain-30): 0.05: max(srain-30); fori=l: length(bin) g(i)=sum(srain-30<=bin(i))./length(srain-30); end x=parmhat(l)*(bin./parmhat (2)); y=-parmhat(l).*log(l-g); plo^Ky/ko*) holdon; plot(0.5: 0.01: 1.5,log(l+(0.5: 0.01: 1.5))); xlabel(\xi(yAsigma)'); ylabel(*-\xilog(l・H(y))'); 结果 Sparmci Sparmhat [-0.0139,5.7800: 0.3829,9.5774] [0.1845,7.4403] 该结果与coles的结果类似. 还可利用evim包进行极值统计,evim包的下载地址: http: //www・sfu.ca/^rgencay/evinLhtml•evim包涉及到这_方面的函数为: out=gpd(datazthreshold,nextremesz9information9); wheredataistheinputvectorandthresholdistheexceedancethresholdvalue.Theoptionnextremessetsthenumberofextremes・Noticethateithernextremesorthresholdshouldbeprovided.Theoptioninformationisastringinputthatdetermineswhetherstandarderrorswillbecalcuhtedwithdataorwillbepostulatedvalue.Itcanbesettoeitherexpectedorobserved,thedefaultvalueisobserved・Outputoutisastructurewiththefollowingsevenfields •uuL.paxeaLb: vcctuiofc^tuiiutrdpuiunirln^thrfkslrlrmrutasc^tiuutedfandlhesrcundisestimatedB •out•funval: valueofthenegativelog-likelihoodfunction. •out•terminated: terminatioacondition.: 1ifsuccessfullyterminated,0otherwise •out.details: detailsofthenonlinearmirumizationprocessofthenegativelikelihood •out•varcov: estimatedvanance-covariancematrixoftheparameters •out•parses: estimatedstandarddeviationsoftheparameters •out.data: vectorofexceedances 程序 clear;clc; load(portpirie.txt1); out=gpd(rain/30/[]/observed>); %outpar_ests(l): shapepara. %out.par_ests (2): scalepara. %95%confidenceintervals: ci_plus=outpar_ests+l・96・\xitpa「_ses; ci_neg=outpa『_ests-L96.★out.parses; %outparses: stdofpara・ 结果 田par_estsfflfurival 田terminated ®details 田varcov 田par_ses 田threshold 田data 田p」ess_threshSci_neg±)ci_plus 2、回归水平的估计当"0时, [0.1845,7.4403] 485.0937 -2 <1x1struct> [0.0102,-0.0655-0.0655.0.9188] [0.1012,0.9585] 30 <152x1double> 0.9913 [-0.0139,5.5615] [0.3829,9.3190] 其中,P{X>zp}=p(“>“)为方便起见,上式中的参数$和&均去掉了头上的帽子, 下文类同,且: gu=P{X>u}^-,其中,R表示“个观测数据中大于“的个数。 n 当Q0时, Z=u+alog— 5 画乙关于log丄的图。 显然,在这种图中,町以较为方便的看出形状参数的符号。 P F图中的直线■根据方框中的公式计算的。 也可以直接根据样本来画returnlevel图。 将犬于“的数据从大到小排列,即: 知S知「••◎,如所对应的回归间隔为: 11+1 注意匕式中”是数据的总数目,N只是人于“的数据数目。 如-T的图如卜图圆点所示。 卜图屮的虎线是95%信度区间,计算方法为: E-1•96严g),°+1•96严口). 其中: %;=[6卅屮-b严{(加久尸-1}+b「(眄jiog(^M)r1{g产-1}]‘ 其中,〃心丄,注盘到此时方差•协方差矩阵的对角线顺序为: P 乙1一「)/刃00' V=0Eyg) 0E@q)Var(a)^ 右下角部分是参数的方差和协方差矩阵。 卜图是以天为单位的returnlevel图: 卜图是以年为单位的returnlevel图: 50 88 returnlevel 223 888 8S88 -00 o52SE5penod(yaar) s functionr$.uml-eve_lgpd(parNCO %reiuren-eve-p-oi %parofGPDfittedi。 data %par(l=shapepara. %par(2=sca-epara・ %>vcosvariance-covariancematrix %dasrQorigma-dasr dasrHdaaro(dasrovihresho-d)j NlLength(da£r)j %empirica-rsurn-eve-empIr-Hsosdasr)jemplrpHl・、((NFl」)・、(-ength(daaro)+l)r %modelreturnlevelx=0.0000005: 0.000001: 0.00& ku=sum(dataO>threshold)./length(dataO); mod_rl=threshold+(par (2)./par(l)).*((ku./x).Apar(l)-l); mod_rp=l./x; %varianceofreturnlevel. new_vcov=zeros(33); new_vcov(l,l)=(ku*(l-ku))/lefigth(dataO); new_vcov(232: 3)=vcov; fori=l: length(x) m=l./xG); dzp=[par⑵*m.Apar(l)*ku.A(par⑴-l);-par (2)*par(l).A(-2)*((m*ku).Apar (1)-1)+p ar⑵*par(l).A(-1)*(m.*ku).Apar(l): *log(m.*ku);par(l).八(W(m.*ku).人par⑴-1)];varzp(i)=dzp'*new_vcov*dzp; clearypdzp; end %95%confidenceinterval posi-ci=mod_rl+1.96*sqrt(varzp); neg_ci=mod」l・L96J*sqrt(varzp); figure(l) semilogx(emp_rp£mp_rI,k.J; holdon; plot(mod_rp,mod_rl); plot(mod_rp,posi_ci/k—7LineWidthrl); plot(mod_rpFneg-ci/k--,/LineWidthrl); xlabelfreturnperiod(day)'); ylabel('returnlevel); figure (2) semilogxCemp^rp./365,emp^rl/k/); holdon; plot(mod_rp./365,mod_rl); plot(mod_rp./365zposi_ci,k-/LineWidth',1); plot(mod_rp./365/neg_ci/k—'/LineWidth',! ); xlabelfreturnperiod(year)1); ylabelfreturnlevel); 3、模型的诊断 经验CDF为: 根据拟合的GPD计算的CDF为: 叽)=1-(1+学)' 画/}(几丿(儿*的图即为P-P图。 程序 functiongpd^probj^lotCpacdataO.threshold)%p-pplot%parofGPDfittedtodata %par(l): shapepara. %par (2): scalepara. %dataO: originaldatadata=dataO(dataOthreshold); N=length(data); %estimateofprobability p2=gpcdf(sort(data-threshold)rpar(l)rpar (2)); plot(plrp2/k.');holdon; plot(0: lr0: l) xlabelfEmpirical); ylabelfModel) QQ图 样本quantde即为x⑴=划+",其对应的累枳概率为: "(N+1),代入GPD分布, 可得模型的quantile为: 画九厂丘“的图即为QQ图。 注意这里的quantile与前面的回归水平爲是不同的,前者 是分布P{X>x\X>“}的quantile,后者是P{X>z}的quantile 程序 functiongp^quan^plotCpaGdataO.threshold) %q-qplot%parofGPDfittedtodata %par(l): shapepara・ %par (2): scalepara. %dataO: originaldata data=dataO(dataOthreshold); N=length(data); %quantileofmodel x=(N: -l: l)./(N+l); q=threshold+(par (2)./par(l));*(x•人(・pai■⑴)-1); plotCq,sort(data),k.');holdon; plot(min(data): 0.01: max(data)rmin(data): 0.01: max(data)) xlabelfModel); ylabelfEmpirical)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 利用 Matlab 进行 极值 统计 一个 例子 GPD 方法