count-dataextreme valuepoisson distribution
I am aware of extreme value theory for continuous distributions. I need to fit an extreme value distribution to the maximum observation of number of events on a day, per month. This seems to be the block maxima problem, which is approximated by the GEV family of distributions for continuous distributions. How do I do this for count data?
As a secondary question, let's assume the basic count process is ~ Poisson. Then does this lead to a different answer to the original question?
To discuss about the constants $a_n$ and $b_n$ in the Fisher-Tippett-Gnedenko theorem, let us assume that $X_k$ is a sequence of i.i.d. r.vs with distribution $F(x)$; we consider $M_n:= \max\{X_1,\,X_2,\,\dots,\,X_n\}$ so that $F_{M_n}(x) = F(x)^n$. The constants are such that $(M_n - b_n) / a_n$ converges to a non-degenerate distribution, say $G(x)$, or in other words: $F(a_n \,x + b_n)^n \to G(x)$ for all $x$. Note that the constants are not uniquely determined: sequences $a_n'$ and $b_n'$ with $a_n' \sim a_n$ and $(b_n - b_n') / a_n \to 0$ can be used as well, with an unchanged non-degenerate limit distribution.
As a general rule, one makes use of the the so-called tail-quantile function $U(t)$. When $F(x)$ is continuous, $U(t)$ is defined for $t > 1$ with its value given by $$ 1 - F(U) = 1 / t. $$ In the vocabulary of applications, $U(t)$ is nothing but the $t$-years return level; the scalar $U$ then has the same physical dimension as the r.vs $X_i$ (length, time, temperature, ...). Now $U(n)$ gives an order of magnitude of $M_n$. Indeed with $U_n:= U(n)$ $$ F_{M_n}(U_n) = F(U_n)^n = \left\{ 1 - [1 - F(U_n)] \right\}^n = \left\{ 1 - 1 / n \right\}^n \approx e^{-1}, $$ so clearly $U(n)$ is in the bulk of the distribution of $M_n$ for large $n$. By either subtracting $U(n)$ or by dividing by $U(n)$ we can hope to scale $M_n$ so that it nearly has a fixed distribution for large $n$. But the choice depends on the type of tail i.e. of the domain of attraction of $X$.
The two simple cases are Weibull and Fréchet. Indeed with $\omega$ being the upper end-point of $F$, it can be proved that \begin{equation} \text{Weibull, type III} \qquad \frac{M_n - \omega}{U(n) - \omega} \to G(x), \end{equation} and \begin{equation} \text{Fréchet, type II} \qquad \frac{M_n - 0}{U(n)} \to G(x). \end{equation} The Gumbel case is more complicated and quite subtle. The value $U(n)$ is then subtracted to $M_n$, i.e. used as $b_n$ and we need the scale $a_n$. If $X$ turns out to have a density $f(x)$ at least near $\omega$, we can use the hazard rate $h(x)$ and the mean residual life $e(x)$ $$ h(x) := \frac{f(x)}{1 - F(x)}, \quad e(x) := \frac{\int_x^\omega [1 - F(t)] \, \text{d}t}{1 - F(x)}. $$ Both $1/h(x)$ and $e(x)$ have the same physical dimension as $x$. Under some mild restrictions we have then \begin{equation} \label{eq:Gum} \text{Gumbel, type I} \qquad \frac{M_n - U(n)}{e(U(n))} \to G(x) \end{equation} One condition ensuring that this convergence holds is one of Von Mises' conditions: the derivative of $1 / h(x)$ exists and tends to $0$ for $x \to \omega$. For many application cases, $h(x)$ is positive and monotonic for $x$ close enough to $\omega$. Assuming this, it can be shown that $F$ is in the Gumbel domain of attraction if and only if the product $h(x)\times e(x)$ tends to $1$ when $x \to \omega$, and then $e(U_n)$ can be replaced by $1 / h(U_n)$.
For many classical distributions such as normal or gamma, neither $F(x)$ nor $U(t)$ are available in closed form. Equivalent quantities can be used, but this requires some math.
As a final remark, note that we can have a finite upper end-point $\omega < \infty$ and yet $F$ in the Gumbel domain of attraction. An example is provided by the reversed Fréchet distribution, for which the determination of the constants is a good exercise.
Disclaimer: At points in the following, this GROSSLY presumes that your data is normally distributed. If you are actually engineering anything then talk to a strong stats professional and let that person sign on the line saying what the level will be. Talk to five of them, or 25 of them. This answer is meant for a civil engineering student asking "why" not for an engineering professional asking "how".
I think the question behind the question is "what is the extreme value distribution?". Yes it is some algebra - symbols. So what? right?
Lets think about 1000 year floods. They are big.
When they happen, they are going to kill a lot of people. Lots of bridges are going down.
You know what bridge isn't going down? I do. You don't ... yet.
Question: Which bridge isn't going down in a 1000 year flood?
Answer: The bridge designed to withstand it.
The data you need to do it your way:
So lets say you have 200 years of daily water data. Is the 1000 year flood in there? Not remotely. You have a sample of one tail of the distribution. You don't have the population. If you knew all of the history of floods then you would have the total population of data. Lets think about this. How many years of data do you need to have, how many samples, in order to have at least one value whose likelihood is 1 in 1000? In a perfect world, you would need at least 1000 samples. The real world is messy, so you need more. You start getting 50/50 odds at about 4000 samples. You start getting guaranteed to have more than 1 at around 20,000 samples. Sample doesn't mean "water one second vs. the next" but a measure for each unique source of variation - like year-to-year variation. One measure over one year, along with another measure over another year constitute two samples. If you don't have 4,000 years of good data then you likely don't have an example 1000 year flood in the data. The good thing is - you don't need that much data to get a good result.
Here is how to get better results with less data:
If you look at the annual maxima, you can fit the "extreme value distribution" to the 200 values of year-max-levels and you will have the distribution that contains the 1000 year flood-level. It will be the algebra, not the actual "how big is it". You can use the equation to determine how big the 1000 year flood will be. Then, given that volume of water - you can build your bridge to resist it. Don't shoot for the exact value, shoot for bigger, otherwise you are designing it to fail on the 1000 year flood. If you are bold, then you can use resampling to figure out how much beyond on the exact 1000 year value you need to build it to in order to have it resist.
Here is why EV/GEV are the relevant analytic forms:
The generalized extreme value distribution is about how much the max varies. The variation in the maximum behaves really different than variation in the mean. The normal distribution, via the central limit theorem, describes a lot of "central tendencies".
Procedure:
now plot the distribution of the result
#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")
This is NOT the "standard normal distribution":
The peak is at 3.2 but the max goes up toward 5.0. It has skew. It doesn't get below about 2.5. If you had actual data (the standard normal) and you just pick the tail, then you are uniformly randomly picking something along this curve. If you get lucky then you are toward the center and not the lower tail. Engineering is about the opposite of luck - it is about achieving consistently the desired results every time. "Random numbers are far too important to leave to chance" (see footnote), especially for an engineer. The analytic function family that best fits this data - the extreme value family of distributions.
Sample fit:
Let's say we have 200 random values of the year-maximum from the standard normal distribution, and we are going to pretend they are our 200 year history of maximum water levels (whatever that means). To get the distribution we would do the following:
(code presumes the above have been run first)
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)
This gives results:
> gevmle(y2)
loc scale shape
3.0965530 0.2957722 -0.1139021
These can be plugged into the generating function to create 20,000 samples
y3 <- rgev(20000,loc=myfit[1],scale=myfit[2],shape=myfit[3])
Building to the following will give 50/50 odds of failing on any year:
mean(y3)
3.23681
Here is the code to determine what the 1000 year "flood" level is:
p1000 <- qgev(1-(1/1000),loc=myfit[1],scale=myfit[2],shape=myfit[3])
p1000
Building to this following should give you 50/50 odds of failing on the 1000 year flood.
p1000
4.510931
To determine the 95% upper CI I used the following code:
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)
The result was:
> mytarget
95%
4.812148
This means, that in order to resist the large majority of 1000 year floods, given that your data is immaculately normal (not likely), you must build for the ...
> out <- pgev(4.812148,loc=fit_d[1],scale=fit_d[2],shape=fit_d[3])
> 1/(1-out)
or the
> 1/(1-out)
shape
1077.829
... 1078 year flood.
Bottom lines:
Best of luck
PS:
PS: more fun - a youtube video (not mine)
https://www.youtube.com/watch?v=EACkiMRT0pc
Footnote: Coveyou, Robert R. "Random number generation is too important to be left to chance." Applied Probability and Monte Carlo Methods and modern aspects of dynamics. Studies in applied mathematics 3 (1969): 70-111.
Best Answer
For any recent visitors, there's been new developments in this area by Hitz, Davis and Samorodnitsky (arXiv:1707.05033). Taking a peaks-over-threshold approach instead of block maxima, the Discrete Generalised Pareto Distribution is derived as the $\operatorname{floor}$ of a GPD, and discrete Maximum Domains of Attraction (DMDA) are introduced by relating them to the classical MDAs. The whole thing is linked to, but different from, Zipf's Law.
In terms of the paper's terminology, the Poisson distribution is in the DMDA of a Gumbel distribution $(\xi = 0)$, as are the Negative Binomial and Geometric distributions.