Extremes

Lesson

Introduction
Theory
Extreme value theory
Multifractal phase transition
To go further
Exercises
- Quiz: GEV distribution
- Quiz: First order phase transition
Conclusion

GEV distribution

The goal of this exercise is to learn who to generate a random variable according to some probability law and plot the corresponding PDF.

Question

Generate a random variable according to the Gumbel distribution with location and scale

Plot the PDF of random variable and the corresponding histogram.

Change the parameters and comment the obtained result.

Do the same thing for Weibull variable

Solution

The Gumbel random variable (RV) can be generated using the Python's function random.gumbel with the appropriate parameters.

In order to plot the PDF, first we calculate the histogram of the RV, then the PDF in computed according to:

$\text{\textbf{Gumbel:}} \\ \quad\\ f(x) = \frac{1}{\sigma} exp(- z - exp(-z)), \quad z = \frac{x- \mu}{\sigma}\\ \quad\\ \text{\textbf{Weibull:}} \\ \quad\\ f(x) = \frac{\xi}{\sigma} \Big( \frac{x}{\sigma} \Big)^{\xi-1} exp(- x/ \sigma)^{\xi}$

"""
@author: yacine.mezemate
"""
import numpy as np
import matplotlib.pylab as plt
#Gumbel
mu, sigma = 0, 0.1 # location and scale
s = np.random.gumbel(mu, sigma, 1000)
plt.figure(0)
count, bins, ignored = plt.hist(s, 30, normed=True)
plt.plot(bins, (1/sigma)*np.exp(-(bins - mu)/sigma) * np.exp( -np.exp( -(bins - mu) 
/sigma) ),lw=3, color='r')
plt.xlabel('Observations')
plt.ylabel('Gumbel PDF')
#Weibull
a = 5. # shape
s = np.random.weibull(a, 1000)
x = np.arange(1,100.)/50.
def weib(x,n,a):
    return (a / n) * (x / n)**(a - 1) * np.exp(-(x / n)**a)
plt.figure(1)
count, bins, ignored = plt.hist(np.random.weibull(5.,1000))
scale = count.max()/weib(x, 1., 5.).max()
plt.figure(1)
plt.plot(x, weib(x, 1., 5.)*scale, lw =3)
plt.xlabel('Observations')
plt.ylabel('Weibull PDF')

"""

@author: yacine.mezemate
"""

import numpy as np
import matplotlib.pylab as plt


#Gumbel
mu, sigma = 0, 0.1 # location and scale
s = np.random.gumbel(mu, sigma, 1000)
plt.figure(0)
count, bins, ignored = plt.hist(s, 30, normed=True)
plt.plot(bins, (1/sigma)*np.exp(-(bins - mu)/sigma) * np.exp( -np.exp( -(bins - mu) 
/sigma) ),lw=3, color='r')
plt.xlabel('Observations')
plt.ylabel('Gumbel PDF')


#Weibull
a = 5. # shape
s = np.random.weibull(a, 1000)
x = np.arange(1,100.)/50.

def weib(x,n,a):
    return (a / n) * (x / n)**(a - 1) * np.exp(-(x / n)**a)

plt.figure(1)
count, bins, ignored = plt.hist(np.random.weibull(5.,1000))
scale = count.max()/weib(x, 1., 5.).max()

plt.figure(1)
plt.plot(x, weib(x, 1., 5.)*scale, lw =3)
plt.xlabel('Observations')
plt.ylabel('Weibull PDF')

(a) : Gumbel distribution; (b): Weibull distribution

Homepage