其中.y>0,l+0,并且广义Pawt。
分布和广义极值分布的参数关系为:
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）
升级会员 