极值分析介绍（2）——分布拟合

发表于 分类于 阅读次数： Disqus：

在第一部分，我们用GEV分布对北京的年最高气温进行了拟合。在这一部分，我们将采用更多的分布进行拟合，并采用KS检验选择其中最优的一个。

读取数据

1
2
3
4
5
6
7
8
9
10
11
12
import pandas as pd
import numpy as np
from plotnine.ggplot import *
from plotnine.qplot import *
from plotnine.geoms import *
from plotnine.coords import *
from plotnine.labels import *
from plotnine.facets import *
from plotnine.scales import *
from plotnine.themes import *

%matplotlib inline
1
2
3
4
5
6
data = pd.read_csv('CHM00054511.csv')
print(data.head())
print(data.dtypes)
data.describe()


   YYYY  MM  DD  TMAX  TMIN  PRCP    SNOW    SNWD
0  1951   1   1  -3.8 -14.0   0.0 -9999.0 -9999.0
1  1951   1   2  -1.2 -10.0   0.5 -9999.0 -9999.0
2  1951   1   3  -4.1 -11.4   0.0 -9999.0 -9999.0
3  1951   1   4  -4.9 -10.1   1.2 -9999.0 -9999.0
4  1951   1   5  -6.6 -10.4   4.5 -9999.0 -9999.0
YYYY      int64
MM        int64
DD        int64
TMAX    float64
TMIN    float64
PRCP    float64
SNOW    float64
SNWD    float64
dtype: object

极值分析

1
2
import lmoments3
from lmoments3 import distr
1
2
3
4
5
6
7
df = data.groupby('YYYY', as_index=False).TMAX.max()
df = df[df.TMAX > 30]
# qplot(x=df.YYYY, y=df.TMAX)
(ggplot(df) +
 aes('YYYY', 'TMAX') +
 geom_point()
)

output_6_0

1
lmoments3.lmom_ratios(df.TMAX, nmom=5)
[37.27076923076924,
 1.0495192307692305,
 0.044338285016250305,
 0.11816135435599627,
 0.004419288498198655]

1
2
3
4
5
6
7
8
gevfit = distr.gev.lmom_fit(df.TMAX)
expfit = distr.exp.lmom_fit(df.TMAX)
gumfit = distr.gum.lmom_fit(df.TMAX)
weifit = distr.wei.lmom_fit(df.TMAX)
gpafit = distr.gpa.lmom_fit(df.TMAX)
pe3fit = distr.pe3.lmom_fit(df.TMAX)
gamfit = distr.gam.lmom_fit(df.TMAX)
glofit = distr.glo.lmom_fit(df.TMAX)

得到不同分布下，在各个回归期的极端最高气温

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
# return years (1.1 to 1000)
T = np.arange(0.1, 999.1, 0.1) + 1

# 估计极端最高气温
fitted_gev = distr.gev(**gevfit)
gevST = fitted_gev.ppf(1.0-1./T)

fitted_exp = distr.exp(**expfit)
expST = fitted_exp.ppf(1.0-1./T)

fitted_gum = distr.gum(**gumfit)
gumST = fitted_gum.ppf(1.0-1./T)

fitted_wei = distr.wei(**weifit)
weiST = fitted_wei.ppf(1.0-1./T)

fitted_gpa = distr.gpa(**gpafit)
gpaST = fitted_gpa.ppf(1.0-1./T)

fitted_pe3 = distr.pe3(**pe3fit)
pe3ST = fitted_pe3.ppf(1.0-1./T)

fitted_gam = distr.gam(**gamfit)
gamST = fitted_gam.ppf(1.0-1./T)

fitted_glo = distr.glo(**glofit)
gloST = fitted_glo.ppf(1.0-1./T)

与经验分布比较

1
2
3
4
5
6
7
8
9
10
11
12
13
14
fitted_dict = {
    'T': T,
    'gevST': gevST,
    'gumST': gumST,
    'expST': expST,
    'weiST': weiST,
    'gpaST': gpaST,
    'pe3ST': pe3ST,
    'gamST': gamST,
    'gloST': gloST
}
fitted_df = pd.DataFrame(fitted_dict)
fitted_df = pd.melt(fitted_df,id_vars='T')
fitted_df.head()
T variable value
0 1.1 gevST 34.854597
1 1.2 gevST 35.454678
2 1.3 gevST 35.847042
3 1.4 gevST 36.144057
4 1.5 gevST 36.384550 
1
2
3
4
5
p = (ggplot(fitted_df) +
 geom_line( aes(x='T', y='value', color='variable')) +
 scale_x_log10()
)
p

output_13_0

1
2
3
4
5
6
7
8
9
10
11
12
N = np.r_[1:len(df.index)+1]*1.0 #must *1.0 to convert int to float
Nmax = max(N)
y = sorted(df.TMAX, reverse=True)
empirical = pd.DataFrame({
    'T':Nmax/N,
    'empirical':y
})
empirical = pd.melt(empirical,id_vars='T')
p = (p + geom_point(aes(x='T', y='value',shape='variable'), empirical,color='orange',show_legend=True) +
scale_color_discrete(name='fitted distribution') +scale_shape_discrete(name='empirical distribution'))
%matplotlib qt5
p

output_14_1

感慨：plotnine作图虽然漂亮，但是在操纵性上还好是matplotlib更顺手一些。

KS检验选择最优的分布

KS检验简介

通常，采用假设检验的方法选择最优分布，在这里我们使用柯尔莫哥洛夫-斯米尔诺夫检验（Колмогоров-Смирнов检验，kolmogorov-smirnov test(K-S test)。KS检验是基于累积概率分布函数，用于检验两个检验分布是否不同或一个经验分布与另一个理想分布是否不同，在Python中可以在scipy.stats.ks_2samp (https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ks_2samp.html)中找到。

如上所述，KS统计量用来量化样本的经验分布函数与参考分布的累积分布函数之间或两个样本的经验分布函数之间的距离。该统计量的零分布是在零假设下计算的，即假设从参考分布中抽取样本（在单样本情况下）或样本是从同一分布中抽取的（在双样本情况下）。在每种情况下，在零假设下考虑的分布是连续分布，但在其他方面不受限制。

双样本K-S检验是用于比较两个样本的最有用和一般的非参数方法之一，因为它对两个样本的经验累积分布函数的位置和形状的差异敏感。KS检验也可以用来做来拟合优度检验（serve as a goodness of fit test）。

双样本K-S检验

假定有分别来自两个独立总体的两个样本，想要检验它们总体分布相同（即零假设）。

  • 参数X，Y，两个由样本观测值组成的1维数据向量，数据的长度可以不同；
  • 输出K-S统计量或者p-value，如果K-S统计量很小或者p-value很大，那么不能拒绝原假设。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
from scipy import stats

# prepare probabilities accroding to observations
N    = np.r_[1:len(df.index)+1]*1.0 #must *1.0 to convert int to float
Nmax = max(N)
P0 = (N-1)/Nmax
P = np.delete(P0,0)
obs = sorted(df.TMAX)

# pextreme high temperatures

gevSTo = fitted_gev.ppf(P)

expSTo = fitted_exp.ppf(P)

gumSTo = fitted_gum.ppf(P)

weiSTo = fitted_wei.ppf(P)

gpaSTo = fitted_gpa.ppf(P)

pe3STo = fitted_pe3.ppf(P)

gamSTo = fitted_gam.ppf(P)

gloSTo = fitted_glo.ppf(P)

# do ks test
ks = [('GEV', stats.ks_2samp(obs, gevSTo)), ('EXP', stats.ks_2samp(obs, expSTo)),
      ('GUM', stats.ks_2samp(obs, gumSTo)), ('WEI', stats.ks_2samp(obs, weiSTo)),
      ('GPA', stats.ks_2samp(obs, gpaSTo)), ('PE3', stats.ks_2samp(obs, pe3STo)),
      ('GAM', stats.ks_2samp(obs, gamSTo)), ('GLO', stats.ks_2samp(obs, gloSTo))]

labels = ['Distribution', 'KS (statistics, pvalue)']
pd.DataFrame(ks, columns=labels)
Distribution KS (statistics, pvalue)
0 GEV (0.07211538461538464, 0.9945953083281929)
1 EXP (0.1942307692307692, 0.15545877513797815)
2 GUM (0.1317307692307692, 0.5996015052982722)
3 WEI (0.07211538461538464, 0.9945953083281929)
4 GPA (0.09639423076923076, 0.9115269503296325)
5 PE3 (0.07283653846153848, 0.9938564597747104)
6 GAM (0.08774038461538464, 0.9568412059042808)
7 GLO (0.08846153846153848, 0.9537725340776948)

可以看到，除了EXP和GUM外，其余的分布的值都在0.9以上，其中GEV（广义极值分布）具有最好的拟合优度（最小的ks统计量和最大的p值）。基于广义极值的结果重新作图(这次不想用pltonine了）:

1
2
3
4
5
6
7
8
a = pd.DataFrame({
'T':1/P,
'value':gevSTo[::-1]})
(ggplot() +
 geom_point(aes(x='T', y='value'),color='orange', data=empirical) +
 geom_line(aes(x='T', y='value'), data=a) +
 scale_x_log10(limits=[1, 100])
)

output_20_0

1
2
3
4
import matplotlib.pyplot as plt
plt.xscale('log')
plt.plot(1/P, gevSTo[::-1])
plt.scatter(Nmax/N, empirical.value, color = 'orangered', facecolors='none', label='Empirical')
<matplotlib.collections.PathCollection at 0x1c11b0ffd0>

output_22_1

1
gevfit
OrderedDict([('c', 0.20571249941146308),
             ('loc', 36.55285724490562),
             ('scale', 1.772333806935174)])

小结

在前一节工作的基础上，我们用多个分布进行了极值分析，并且利用KS检验选择了其中最好的一个。接下来，我们将采用Bootstrap采样方法估计选择的分布的置信区间。