极值分析介绍（1）——入门

发表于 2019-01-28 分类于他山攻错阅读次数： Disqus：

以前做多元统计的时候，主要关注正常数据而非极端数据的分布，极值分析刚好相反，算是给我一个新的视角，所以做了一些学习。

介绍

极值分析（Extreme value analysis, EVA）,是一个主要用来估计比历史观测到的事件更加极端的事件的概率的方法，广泛用于金融、工程、气象、水文、海洋等诸多领域。

EVA就不用我多说了吧…

基于Python对时间序列（time series）数据进行分析包括以下基本步骤。

准备数据（Prepare Data）

中国气象数据网

以气象数据为例，我国的气象数据可以通过中国气象数据网（http://data.cma.cn ）获取。

中国气象数据网是气象科学数据共享网的升级系统，是国家科技基础条件平台的重要组成部分，是气象云的主要门户应用系统，是中国气象局面向国内和全球用户开放气象数据资源的权威的、统一的共享服务平台，是开放我国气象服务市场、促进气象信息资源共享和高效应用、构建新型气象服务体系的数据支撑平台。

例如，较为常用的"中国地面气候资料日值数据集（V3.0）"包含了中国699个基准、基本气象站1951年1月以来本站气压、气温、降水量、蒸发量、相对湿度、风向风速、日照时数和0cm地温要素的日值数据。但是，要求教育科研实名注册用户才能访问。

data-cma

基于GHCNpy下载NOAA数据

山重水复疑无路，柳暗花明又一村。可以从NOAA的网站（https://www.ncdc.noaa.gov ）下载数据，具体而言是，全球历史气候学网络——逐日数据集（The Global Historical Climatology Network-Daily Data Set, GHCN-D），GHCN-D在日常尺度上为地球气候提供了坚实的基础，并且是美国每日天气数据的官方档案。数据集每晚更新，新数据的滞后时间约为一天。

由于边境墙的问题，美国政府关门了，结果导致NOAA的网站无法访问，好在FTP还可以照常访问。网上有下载GHCN的Python包，我选择了GHCNPy。麻烦的是GHCNpy的最新版本还是3年前的，并且是基于Python2的，而我用的是Python 3，经过一番折腾之后，姑且先凑活着用了。（找到了一个更加活跃也更加新的包，不过先不用了，这不是重点。）

noaa-shutdown

首先，查找目标城市可以获取的站点，以北京为例：

import ghcnpy as gp
import ghcnpy.get_station_id as gsi
import ghcnpy.get_dly as gd

1	gp.find_station('Beijing')

LOOKUP BY STATION NAME:  Beijing

GRABBING LATEST STATION METADATA FILE
GHCND ID          LAT        LON    ELEV  ST       STATION NAME
###############################################################
CHM00054511   39.9330   116.2830    55.0      BEIJING

第二步，根据站点的名称下载数据，并保存为csv格式：

1	gp.output_to_csv("CHM00054511")

OUTPUTTING TO CSV:  CHM00054511 .csv

GETTING DATA FOR STATION:  CHM00054511

利用Pandas处理时间序列

作为数据科学中非常强大的一个库，Pandas提供了一组标准的时间序列处理工具和数据算法，可以高效处理非常大的时间序列，对于金融和经济数据有为有用。言归正传，首先导入需要的库

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 *

1	%matplotlib inline

导入数据，并查看数据的基本情况：

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

取每一年的最大值作为极值

df = data.groupby('YYYY', as_index=False).TMAX.max()
# qplot(x=df.YYYY, y=df.TMAX)
(ggplot(df) +
 aes('YYYY', 'TMAX') +
 geom_point()
)

output_19_0

极值分析示例

L-moments介绍

在统计学中，线性矩（L-moments）是用于概括描述概率分布形状的一系列统计量。他们是与常规矩相似的阶次统计量（L-statistics, L代表linear)的的现行组合，可以用于计算类似于标准差、偏度和峰度的量，分别称为L-尺度（L-scale)、L-偏度（L-skewness）和L-峰度（L-kurtosis）（注意：L-mean与传统的均值相同）。同样类比于标准化矩，标准化线性矩（Standardised L-moments)被称为线性矩比（L-moments ratio)。

正如传统的矩估计一样，一个概率分布具有一组总体线性矩（popluation L-moments)。对于从总体中得到的样本，可以定义样本线性矩（sample L-moments)，并且可以作为总体线性矩的估计量。

以上翻译自维基百科（https://en.wikipedia.org/wiki/L-moment ），如有错误之处，欢迎指出，感激不尽！

Python中计算L-moments的包为lmoments(for Python 2.x)和lmoments3(for Python 3.x)，可以用来拟合分布的参数。lmoments3(https://pypi.org/project/lmoments3/ )支持许多概率分布，具体包括以下(括号内为lmoments3中的名称）：

Exponential(exp)
Gamma(gam)
Generalized Extreme Value(gev)
Generalized Logistic(glo)
Generalized Normal(gno)
Generalized Pareto(gpa)
Gumbel(gum)
Kappa(kap)
Normal(nor)
Pearson III(pe3)
Wakeby(wak)
Weibull(wei)

上述分布可以通过引用loments3.distr模块使用。

以广义极值分布（GEV）为例

导入lmoments包，计算l-moments

1 2	import lmoments3 from lmoments3 import distr

1	lmoments3.lmom_ratios(df.TMAX, nmom=5)

[36.24117647058823,
 1.9677787532923634,
 -0.4057142161813552,
 0.4596155640192579,
 -0.3562413164625048]

1	import lmoments

1	LMU = lmoments.samlmu(df.TMAX)

拟合一个GEV分布

1 2	gevfit = lmoments.pelgev(LMU) print(gevfit)

[36.561393299164436, 3.796087323276009, 1.1988788928571046]

1	paras = distr.gev.lmom_fit(df.TMAX)

选择若干回归期，将回归其转换为概率，得到这些概率所对应的极值

avi  = [2,5,10,20,50,100,200,500,1000];

pavi =[(x - 1) / x for x in avi]

print(pavi)

gevqua = lmoments.quagev(pavi, gevfit)
gevqua

[0.5, 0.8, 0.9, 0.95, 0.98, 0.99, 0.995, 0.998, 0.999]





array([37.68728914, 39.20344246, 39.51451693, 39.63779132, 39.69831662,
       39.71501032, 39.72222143, 39.72591537, 39.72695558])

1 2	fitted_gev = distr.gev(**paras) fitted_gev.ppf(pavi)

<lmoments3.distr.LmomFrozenDistr at 0x26265aa35f8>

所以，根据上面的结果，可以得知，北京2年1遇的极端最高气温为37.687摄氏度，10年1遇的极端最高气温为39.615摄氏度，100年一遇则为39.715摄氏度。

仔细一看，这个结果跟我们的实际观感好像不相符啊，北京这几年动辄上个40度好像也不是太稀奇的事情，难道百年一遇的好事情都被咱们赶上啦？

其实看一下上面的图就发现了，有三年的最高气温还不到30度，最近两年甚至低于10度，明摆着是数据错了，在下一节中，我们将剔除掉异常的数据，可以发现结果会有比较大的差别。

所以，在分析实际数据时，对于数据的理解、数据的质量和清洗对于结果的准确性和合理性非常重要。对于分析的结果也不能盲目相信，还是要有一些经验上的判断。

结果校验

1	qplot(x=avi, y=gevqua, xlab='avi', ylab='gevqua')

output_38_0

小结

上述工作演示了如何对感兴趣的时间序列（北京市历年最高气温）进行初步的极值分析。

实际上，在实践当中，极值分析要复杂得多，其中一个很重要的原因就是高质量的数据并不是随手可得的。我们要花费大量的时间用于数据上，包括数据的规整、缺失值的处理和野值的判断及处理上。

此外，在上述示例中，我们只用了一种分布，而实际当中往往会同时采用多种不同的分布去做拟合,后使用拟合优度度量（例如，Anderson-Darling检验）来选择最佳的分布。

这只是开始，接下来将更加深入地去做些分析。

介绍