Análisis de extremos

Análisis de extremos
Referencias
• Wilks (sección 4.4.5): “dice mucho, explica poco”
• Coles (2001) “An Introduction to Statistical
Modeling of Extreme Values”
Introducción
• Objetivo del análisis de extremos: cuantificar el
comportamiento aleatorio de un determinado
proceso para valores atípicamente altos o bajos.
• 1ª dificultad  se busca estimar la probabilidad de
ocurrencia de eventos mayores (o menores) a los
observados; i.e. se requieren modelos para
“extrapolar” los datos
• Por ejemplo:
Variable aleatoria: “Máxima velocidad de viento
anual”:
Máxima velocidad de viento media en 10 min a 10 m
de altura registrada cada 1 hora en una determinada
estación durante 34 años
25
V [m/s]
20
15
10
0
5
10
15
20
Años
25
30
35
Se puede estimar la distribución empírica como:
25
V [m/s]
20
15
10
0
10
1
10
Tr = 1/(1-F)
2
10
25
V [m/s]
20
15
10
0
10
1
10
Tr = 1/(1-F)
2
10
Período de retorno de un valor Zs es el número esperado
de observaciones de Z entre dos excedencias
consecutivas de Zs.
1
Tr 
1  P( Z  Z S )
Si la serie es de datos anuales entonces Tr está en años.
Interés social: ¿qué valores es esperable que ocurran
con alto período de retorno? (inundaciones, sequías,
olas de calor, vientos…)
Valores objetivo típicos corresponden a entre 100 y
10.000 años.
?
25
V [m/s]
20
15
10
0
10
1
10
Tr = 1/(1-F)
2
10
Introducción
• Objetivo del análisis de extremos: cuantificar el
comportamiento aleatorio de un determinado
proceso para valores atípicamente altos o bajos.
• 2ª dificultad  La información disponible es
escasa; por definición se está trabajando con
valores atípicos (poco probables).
• Por ejemplo:
De una determinada variable tenemos datos horarios
durante varios años
30
20
10
0
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
4
x 10
1500
1000
500
0
0
5
10
15
20
25
Introducción
• 2ª dificultad  La información disponible es
escasa; por definición se está trabajando con
valores atípicos (poco probables).
• Muy importante:
1) Maximizar el uso de la información
2) Calcular intervalos de
Introducción
• Existen dos aproximación “clásicas” al cálculo de extremos:
(1) Análisis de los máximos alcanzados en un período de
tiempo fijo (típicamente máximos anuales)
(2) Análisis de los máximos alcanzados en cada uno de los
“eventos extremos” registrados (picos sobre el umbral)
El análisis (1) considera un “evento extremo” por año de
registro; en el análisis (2) puede haber años con más de un
“evento extremo” y años sin “eventos extremos”.
Máximos anuales
• Coles (2001) capítulo 3 (+ algo de capítulo 2)
Máximos anuales
• Nuestra variable aleatoria va a ser Mn:
con X1,…,Xn una secuencia de variables aleatorias
independientes con la misma distribución F(x).
Típicamente Xi son los valores obtenidos al medir un
determinado proceso con una frecuencia fija (e.g.
datos de viento horarios), y n es la cantidad de datos
en una año (e.g. 365 x 24).
Máximos anuales
• En principio la distribución de Mn podría derivarse
de la distribución F(x):
… pero pequeños errores en el ajuste de la cola
superior se traducen en grandes errores en la
distribución de máximos anuales
• Por ejemplo:
De una determinada variable tenemos datos horarios
durante varios años  se puede ajustar una
distribución F(x)
30
20
10
0
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
4
x 10
1500
20
f(x)
1000
¡?
15
10
5
500
0
18
0
0
5
10
15
20
25
20
22
24
26
28
30
32
Máximos anuales
• ¿Alternativa?  Extremal Type Theorem
Dice que…
…el máximo Mn de n observaciones independientes
de Xi tiende a una distribución de probabilidad G(z)
fija cuando n tiende a infinito, independientemente
de cuál sea la distribución F(x)
La distribución G(z) se denomina distribución
generalizada de extremos (GEV distribution).
Máximos anuales
• En la realidad:
• n no tiende a infinito
• Las observaciones no son independientes
• La distribución F(x) no necesariamente es estacionaria
• A pesar de lo cual el uso de la distribución GEV a probado
ser una herramienta valiosa en el análisis de extremos …
• … siempre y cuando se tenga una comprensión suficiente
de los procesos físicos que subyacen a los datos analizados,
de modo de aplicar la GEV “con criterio”
Máximos anuales
• Distribución GEV
con
siendo:
mu  parámetro de posición
sigma  parámetro de escala
xi
 parámetro de forma
Máximos anuales
• El parámetro de forma (xi) determina si la distribución está acotada o
no…
xi = 0  Distribución Tipo I (Gumbel) no acotada
xi > 0  Distribución Tipo II (Frechet o Fisher-Tippett II) con límite inferior
xi < 0  Distribución Tipo III (Weibull) con límite superior
Máximos anuales
• Cálculo de cuantiles
Parea calcular el cuantil de probabilidad de
excedencia p (zp), se invierte la distribución G(z),
obteniéndose:
Cuando se trabaja con máximos anuales el período
de retorno es Tr = 1/p
Máximos anuales
• Papel de Gumbel:
Si se define la variable
, siendo p la
probabilidad de excedencia de zp (G(zp) = 1-p ), entonces el
gráfico log(yp) – zp será lineal en el caso de la distribución de
Gumbel (xi = 0)
Tipo III: Weibull
Tipo I: Gumbel
Tipo II: Frechet
Máximos anuales
• Ejemplo GEV:
W
Velocidad de viento V [m/s]
25
20
15
10
0
10
1
2
10
10
Período de retorno [años]
3
10
Máximos anuales
• Ajuste de los parámetros de la distribución GEV
(1) Máxima verosimilitud
(2) Probability Weighted Moments (PWM;
combinación lineal de momentos de la
distribución)
Máximos anuales
• Ajuste parámetros GEV mediante máxima
verosimilitud
Da resultados adecuados cuando xi>-0,5 (habitual en
la práctica) y cuando el número de datos disponible
es “grande”
Se buscan los parámetros que maximizan el
logaritmo de la función de verosimilitud.
Máximos anuales
En MATLAB el comando es: gevfit
Máximos anuales
• Ajuste parámetros GEV mediante PWM
Da mejor ajuste cuando el número de datos
disponible es “chico”
Equivalente al método de los momentos, pero se
igualan combinaciones lineales de momentos en
lugar de igualar momentos.
Máximos anuales
• Ajuste parámetros GEV mediante PWM
El método no está programado en MATLAB pero es
de fácil programación.
Máximos anuales
W
Velocidad de viento V [m/s]
25
¿?
20
15
data1
ML
PWM
10
0
10
1
2
10
10
Período de retorno [años]
3
10
Máximos anuales
• Cálculo de intervalos de confianza
(1) Método Delta
(2) Perfil de verosimilitud (profile likelihood)
(3) Métodos basados en Bootstrapping paramétrico
(muestreo aleatorio y simulación)
Máximos anuales
• Intervalos de confianza: Método Delta
Asigna a los cuantiles una distribución normal 
intervalos de confianza simétricos en torno al valor
estimado…
… no es realista cuando se extrapola a altos períodos de
retorno: se corre el riesgo de subestimar límite superior.
No está programado en MATLAB; puede programarse con
cierta facilidad siguiendo Coles (2001, Cap.2)
Máximos anuales
• Intervalos de confianza: Perfil de Verosimilitud
Máximos anuales
• Intervalos de confianza: Perfil de Verosimilitud
Máximos anuales
• Intervalos de confianza: Perfil de Verosimilitud
Máximos anuales
• Intervalos de confianza: Perfil de Verosimilitud
Para usar el perfil de verosimilitud para estimar los
intervalos de confianza de los cuantiles de la GEV es
necesario expresar uno de los parámetros de la
distribución en función del valor del cuantil (ver
Coles 2001, Cap 3).
Máximos anuales
• Intervalos de confianza: Perfil de Verosimilitud
Máximos anuales
• Intervalos de confianza: Perfil de Verosimilitud
No está programado en MATLAB
Es posible programarlo siguiendo Coles (2001)
capítulos 2 y 3, pero requiere usar funciones de
optimización (búsqueda de mínimos) para construir
los perfiles de verosimilitud.
Máximos anuales
• Intervalos de confianza: Bootstrapping paramétrico
• Aproximación tradicional; dada GEV ajustada con
datos originales:
(1) se simulan N series de igual duración que la serie
original
(2) a cada serie se le ajusta una GEV y se calcula el
cuantil de interés
(3) se estiman intervalos de confianza a partir de
datos generados
Máximos anuales
• Intervalos de confianza: Bootstrapping paramétrico
• Aproximación alternativa; dada GEV ajustada con datos
originales:
(1) se estima matriz de covarianza de los parámetros
(2) se simulan N conjuntos de parámetros usando la
matriz de covarianza
(3) Se estima el cuantil objetivo con cada conjunto de
parámetros
(4) se estiman intervalos de confianza a partir de datos
generados
Máximos anuales
• Intervalos de confianza: Bootstrapping paramétrico
No están programados en MATLAB pero son de
sencilla programación una vez que se tiene la matriz
de covarianza
Máximos anuales
• Intervalos de confianza: Bootstrapping paramétrico
Matriz de covarianza de los parámetros para ajuste
ML: función gevlike
Matriz de covarianza de los parámetros para ajuste
PWM hay que programarla usando
Máximos anuales
25
W
Velocidad de viento V [m/s]
30
20
data1
ML
PWM
15
10
0
10
1
2
10
10
Período de retorno [años]
3
10
Máximos anuales
• Ejemplo: análisis de extremos de la altura de ola
significante en la costa de Rocha
x 10
Profundidad (mWh)
6
25
6.21
20
5
6.2
10
15
20
22.5
Y (UTM zona 22H)
6.19
15
25
20
6.18
25
15
20
6.17
ADCP
10
22.5
10
15
6.16
10
20
B11
5
15
25
Condición de borde
6.15
1.9
2
2.1
2.2
2.3
X (UTM zona 22H)
2.4
2.5
2.6
x 10
5
0
Máximos anuales
Máximo Hm0 Anual [m]
6
4
2
0
1980
1985
1990
1995
Año
2000
2005
2010
Máximos anuales
Distribución de extremos GEV de Hm0 (Omnidireccional)
Altura de ola significante espectral Hm0 [m]
9
8
7
6
5
4
3
0
10
1
2
10
10
Período de retorno [años]
3
10
Altura de ola significante espectral Hm0 [m]
8
6
4
2
10
0
10
10
Período de retorno [años]
0
2
SSE
10
8
6
4
2
10
Período de retorno [años]
2
10
8
6
4
2
10
0
10
10
Período de retorno [años]
0
2
S
10
8
6
4
2
10
Período de retorno [años]
2
Altura de ola significante espectral Hm0 [m]
E
Altura de ola significante espectral Hm0 [m]
10
Altura de ola significante espectral Hm0 [m]
ENE
Altura de ola significante espectral Hm0 [m]
Altura de ola significante espectral Hm0 [m]
ESE
10
8
6
4
2
10
0
10
10
Período de retorno [años]
0
2
SSW
10
8
6
4
2
10
Período de retorno [años]
2
Altura de ola significante espectral Hm0 [m]
Máximos anuales
SE
10
8
6
4
2
10
0
10
Período de retorno [años]
2
Máximos anuales
Distribución de extremos GEV de Hm0 (Omnidireccional)
Distribución de extremos GEV de Hm0 (Direccional)
11
Altura de ola significante espectral Hm0 [m]
Altura de ola significante espectral Hm0 [m]
9
8
7
6
5
4
3
0
10
1
2
10
10
Período de retorno [años]
3
10
10
9
8
7
6
5
4
3
0
10
1
2
10
10
Período de retorno [años]
3
10