我来自土木工程,我们使用极值理论,例如GEV分布来预测某些事件的值,例如最大风速,即98.5%的风速会低于该值。
我的问题是为什么要使用这样的极值分布?如果我们只使用整体分布并获得98.5% 概率的值,会不会更容易?
我来自土木工程,我们使用极值理论,例如GEV分布来预测某些事件的值,例如最大风速,即98.5%的风速会低于该值。
我的问题是为什么要使用这样的极值分布?如果我们只使用整体分布并获得98.5% 概率的值,会不会更容易?
免责声明:在以下几点,这大致假定您的数据是正态分布的。如果您实际上是在设计任何东西,那么请与强大的统计专家交谈,并让该人在线上签名,说明级别将是多少。与其中 5 人或 25 人交谈。此答案适用于询问“为什么”的土木工程专业学生,而不适用于询问“如何”的工程专业人士。
我认为问题背后的问题是“什么是极值分布?”。是的,它是一些代数 - 符号。所以呢?对?
让我们想想 1000 年的洪水。他们很大。
当它们发生时,它们会杀死很多人。许多桥梁正在倒塌。
你知道哪座桥不会倒塌吗?我愿意。你还没有……还没有。
问题:哪座桥不会在 1000 年的洪水中倒塌?
答:设计的桥梁可以承受它。
您需要按照自己的方式进行操作的数据:
因此,假设您拥有 200 年的每日水数据。那里有1000年的洪水吗?不是远程。您有一个分布尾部的样本。你没有人口。如果您了解所有洪水的历史,那么您将拥有全部数据。让我们考虑一下。您需要多少年的数据,多少个样本才能获得至少一个可能性为千分之一的值?在一个完美的世界中,您至少需要 1000 个样本。现实世界是混乱的,所以你需要更多。您开始在大约 4000 个样本处获得 50/50 的赔率。您开始得到保证在大约 20,000 个样本中拥有超过 1 个样本。样本并不意味着“一秒钟与下一秒钟的水”,而是对每个独特变异来源的衡量——比如逐年变化。一项措施超过一年,与另一年的另一项措施一起构成两个样本。如果您没有 4,000 年的良好数据,那么您可能在数据中没有 1000 年洪水的示例。好消息是 - 你不需要那么多数据来获得好的结果。
以下是如何用更少的数据获得更好的结果:
如果您查看年度最大值,您可以将“极值分布”拟合到 year-max-levels 的 200 个值,您将得到包含 1000 年洪水的分布-等级。它将是代数,而不是实际的“它有多大”。您可以使用该等式来确定 1000 年洪水的规模。然后,鉴于水的体积 - 你可以建造你的桥梁来抵抗它。不要追求精确的值,而要追求更大的值,否则您将其设计为在 1000 年的洪水中失败。如果您是大胆的,那么您可以使用重新采样来确定您需要将其构建到的确切 1000 年值超出多少才能使其抵抗。
这就是为什么 EV/GEV 是相关分析形式的原因:
广义极值分布是关于最大值变化的程度。最大值的变化与平均值的变化完全不同。正态分布通过中心极限定理描述了许多“中心趋势”。
程序:
现在绘制结果的分布
#libraries
library(ggplot2)
#parameters and pre-declarations
nrolls <- 1000
ntimes <- 10000
store <- vector(length=ntimes)
#main loop
for (i in 1:ntimes){
#get samples
y <- rnorm(nrolls,mean=0,sd=1)
#store max
store[i] <- max(y)
}
#plot
ggplot(data=data.frame(store), aes(store)) +
geom_histogram(aes(y = ..density..),
col="red",
fill="green",
alpha = .2) +
geom_density(col=2) +
labs(title="Histogram for Max") +
labs(x="Max", y="Count")
这不是“标准正态分布”:
峰值为 3.2,但最大值上升至 5.0。它有偏斜。它不会低于大约 2.5。如果您有实际数据(标准法线)并且您只选择尾部,那么您将沿着这条曲线均匀随机地选择一些东西。如果你幸运的话,你会朝着中心而不是下尾巴。工程与运气相反——它是关于每次都能始终如一地实现预期的结果。“随机数太重要了,不能任凭运气”(见脚注),尤其是对于工程师而言。最适合该数据的分析函数族 - 分布的极值族。
样本拟合:
假设我们有 200 个来自标准正态分布的年度最大值的随机值,我们将假装它们是我们 200 年的最高水位历史(无论这意味着什么)。要获得分布,我们将执行以下操作:
(代码假定上面已经先运行)
library(SpatialExtremes) #if it isn't here install it, it is the ev library
y2 <- sample(store,size=200,replace=FALSE) #this is our data
myfit <- gevmle(y2)
这给出了结果:
> gevmle(y2)
loc scale shape
3.0965530 0.2957722 -0.1139021
这些可以插入到生成函数中以创建 20,000 个样本
y3 <- rgev(20000,loc=myfit[1],scale=myfit[2],shape=myfit[3])
建立以下条件将使任何一年的失败几率为 50/50:
平均值(y3)
3.23681
这是确定 1000 年“洪水”水平的代码:
p1000 <- qgev(1-(1/1000),loc=myfit[1],scale=myfit[2],shape=myfit[3])
p1000
建立在以下基础上应该会给你 50/50 的失败几率在 1000 年的洪水中失败。
p1000
4.510931
为了确定 95% 的 CI 上限,我使用了以下代码:
myloc <- 3.0965530
myscale <- 0.2957722
myshape <- -0.1139021
N <- 1000
m <- 200
p_1000 <- vector(length=N)
yd <- vector(length=m)
for (i in 1:N){
#generate samples
yd <- rgev(m,loc=myloc,scale=myscale,shape=myshape)
#compute fit
fit_d <- gevmle(yd)
#compute quantile
p_1000[i] <- qgev(1-(1/1000),loc=fit_d[1],scale=fit_d[2],shape=fit_d[3])
}
mytarget <- quantile(p_1000,probs=0.95)
结果是:
> mytarget
95%
4.812148
这意味着,为了抵御 1000 年的大部分洪水,鉴于您的数据非常正常(不太可能),您必须为...
> out <- pgev(4.812148,loc=fit_d[1],scale=fit_d[2],shape=fit_d[3])
> 1/(1-out)
或者
> 1/(1-out)
shape
1077.829
... 1078 年洪水。
底线:
祝你好运
PS:
PS:更有趣 - 一个 youtube 视频(不是我的)
https://www.youtube.com/watch?v=EACkiMRT0pc
脚注:Coveyou, Robert R. “随机数生成太重要了,不能靠运气。” 应用概率和蒙特卡罗方法以及动力学的现代方面。应用数学研究 3(1969):70-111。
如果您只对尾巴感兴趣,那么将数据收集和分析工作集中在尾巴上是有意义的。这样做应该更有效率。我强调了数据收集,因为在提出 EVT 分布的论点时,这方面经常被忽略。事实上,收集相关数据来估计某些领域的整体分布可能是不可行的。我将在下面更详细地解释。
如果您像@EngrStudent 的示例一样查看千分之一的洪水,那么要构建正态分布的主体,您需要大量数据来填充观察结果。您可能需要过去数百年来发生的每一次洪水。
现在停下来想想到底什么是洪水?当我的后院在一场大雨后被洪水淹没时,是洪水吗?可能不是,但是将洪水与非洪水事件划定的界限到底在哪里?这个简单的问题突出了数据收集的问题。您如何确保我们按照相同的标准收集身体上的所有数据几十年甚至几个世纪?要收集洪水分布体的数据几乎是不可能的。
因此,这不仅仅是分析效率的问题,而是数据收集的可行性问题:是对整个分布建模还是仅对尾部建模?
自然,有了尾巴,数据收集就容易多了。如果我们为大洪水定义足够高的阈值,那么我们就有更大的机会以某种方式记录所有或几乎所有此类事件。很难错过一场毁灭性的洪水,如果有任何文明存在,就会保存一些关于该事件的记忆。因此,鉴于数据收集在极端事件上比在可靠性研究等许多领域中的非极端事件上更加稳健,因此构建专门关注尾部的分析工具是有意义的。
您使用极值理论从观察到的数据中进行推断。通常,您拥有的数据不足以为您提供尾部概率的合理估计。以@EngrStudent 的千分之一年事件为例:这对应于找到分布的 99.9% 分位数。但是,如果您只有 200 年的数据,则最多只能计算 99.5% 的经验分位数估计值。
极值理论让您通过对尾部分布的形状做出各种假设来估计 99.9% 的分位数:它是平滑的,它以某种模式衰减,等等。
您可能会认为 99.5% 和 99.9% 之间的差异很小;毕竟只有0.4%。但这是概率上的差异,当你处于尾部时,它可以转化为分位数的巨大差异。这是伽马分布的示意图,随着这些事情的发展,它没有很长的尾巴。蓝线对应 99.5% 分位数,红线对应 99.9% 分位数。虽然在垂直轴上它们之间的差异很小,但在水平轴上的分离却很大。对于真正的长尾分布,这种分离只会变得更大;伽马实际上是一个相当无害的案例。
通常,基础数据(例如,高斯风速)的分布是针对单个样本点的。第 98 个百分位会告诉您,对于任何随机选择的点,该值有 2% 的机会大于第 98 个百分位。
我不是土木工程师,但我想你想知道的不是任何一天的风速超过某个数字的可能性,而是最大可能阵风的分布,比如说,一年的过程。在这种情况下,如果每日阵风最大值呈指数分布,那么您想要的是 365 天内最大阵风的分布……这就是极值分布要解决的问题。