matlab模拟gpd,利用Matlab进行极值统计的一个例子——GPD方法.PDF

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

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

科研菜鸟

/u/sanshiphy

数据和例子均来自于S. Coles, An introduction to statistical modeling of extreme values,

Springer, 2007

一、 数据

数据是south-west England 这个地方的日降雨量 (1914-1962)，原始数据可以从R 语言

包中调出：

0

8

0

6

n

i

a

r 0

4

0

2

0

0 5000 10000 15000

Index

二、 GPD 的极大似然估计

1、超越量的似然估计

广义Pareto 分布为：



 ξy −1/ξ

− < > ≡ =− +

P {X u y | X u } H (y ) 1 1  

 σ 



其中，y >0 ，1 y /  0 ，并且广义Pareto 分布和广义极值分布的参数关系为：

+ξ σ >





ξ ξ ，σ σ =+ξ(u −µ) 。

Matlab 中涉及这一方面的函数主要有：

parmhat = gpfit(X)

returns maximum likelihood estimates of the parameters for the two-parameter generalized Pareto

(GP) distribution given the data in X. parmhat(1) is the tail index (shape) parameter, K and

parmhat(2) is the scale parameter, sigma. gpfit does not fit a threshold (location) parameter.

[parmhat,parmci] = gpfit(X)

returns 95% confidence intervals for the parameter estimates.

程序

clear;clc;

load('rain.txt');

srain=rain(rain>30);

[parmhat,parmci]=gpfit(srain-30);

bin=min(srain-30):0.05:max(srain-30);

for i=1:length(bin)

g(i)=sum(srain-30<=bin(i))./length(srain-30);

end

x=parmhat(1).*(bin./parmhat(2));

y=-parmhat(1).*log(1-g);

plot(x,y,'ko')

hold on;

plot(0.5:0.01:1.5,log(1+(0.5:0.01:1.5)));

xlabel('xi(y/sigma)');

ylabel('-xilog(1-H(y))');

结果

该结果与Coles 的结果类似。

还可利用evim 包