User guide

First of all you should import the package:

import skextremes as ske
import matplotlib.pyplot as plt
%matplotlib inline

Some datasets are included in the package. For example, we will use sea level data from Port Pirie, in Australia.

data = ske.datasets.portpirie()

print(data.description)

Annual Maximum Sea Levels at Port Pirie, South Australia
--------------------------------------------------------

Fields:
 year: numpy.array defining the year for the row data.
 sea_level: numpy.array defining annual maximum sea level recorded at
    Port Pirie, South Australia.

Source:
 -Coles, S. G. (2001). An Introduction to Statistical Modelling of
    Extreme Values. London: Springer.

As can be seen from the description of the dataset we have two fields, sea_level, with the annual maximum sea level records, and year, indicating the year of the record.

To get the dataset we can use the .asarray method to obtain all the fields as a numpy.array or just select a field to get a 1D numpy.array with the records for the field.

data_array = data.asarray()
sea_levels = data.fields.sea_level

print(sea_levels)

[ 4.03  3.83  3.65  3.88  4.01  4.08  4.18  3.8   4.36  3.96  3.98  4.69
  3.85  3.96  3.85  3.93  3.75  3.63  3.57  4.25  3.97  4.05  4.24  4.22
  3.73  4.37  4.06  3.71  3.96  4.06  4.55  3.79  3.89  4.11  3.85  3.86
  3.86  4.21  4.01  4.11  4.24  3.96  4.21  3.74  3.85  3.88  3.66  4.11
  3.71  4.18  3.9   3.78  3.91  3.72  4.    3.66  3.62  4.33  4.55  3.75
  4.08  3.9   3.88  3.94  4.33]

We have several options, located in the skextremes.models subpackage, to calculate extreme values. Depending the chosen model some methods and attributes will be available.

Models are divided in several packages:

  • wind: Some very basic functions to calculate wind extreme values obtained from the Wind Energy industry. These are very basic approximations based on parameters of wind data and should only be used to calculate the expected wind speed value for a 50 years return period.

  • engineering: Basic models found in the literature based, mainly, on the Gumbel distribution and the Block Maxima approach.

  • classic: A more classical approach from the theory to obtain extreme values using fitting accepted distributions generally used in the extreme value theory.

model = ske.models.engineering.Lieblein(sea_levels)

Once the model is fitted you can see the results:

model.plot_summary()

(<matplotlib.figure.Figure at 0xf449d70>,
 <matplotlib.axes._subplots.AxesSubplot at 0x4632f70>,
 <matplotlib.axes._subplots.AxesSubplot at 0x466a410>,
 <matplotlib.axes._subplots.AxesSubplot at 0x4693410>)
_images/User_guide_12_1.png

To see, for instance, the parameters obtained you can use:

print(model.c, model.loc, model.scale)

0 3.86774852633 0.198420194778

Another approach would be to use a more classic model:

model = ske.models.classic.GEV(sea_levels, fit_method = 'mle', ci = 0.05,
                              ci_method = 'delta')

d:usersX003621AppDataLocalContinuumMiniconda3libsite-packagesnumdifftoolscore.py:753: UserWarning: The stepsize (3.16814) is possibly too large!
  warnings.warn('The stepsize (%g) is possibly too large!' % h1[i])
d:usersX003621AppDataLocalContinuumMiniconda3libsite-packagesnumdifftoolscore.py:753: UserWarning: The stepsize (0.18069) is possibly too large!
  warnings.warn('The stepsize (%g) is possibly too large!' % h1[i])
d:usersX003621AppDataLocalContinuumMiniconda3libsite-packagesnumdifftoolscore.py:753: UserWarning: The stepsize (0.278302) is possibly too large!
  warnings.warn('The stepsize (%g) is possibly too large!' % h1[i])
d:usersX003621AppDataLocalContinuumMiniconda3libsite-packagesnumdifftoolscore.py:753: UserWarning: The stepsize (0.0664633) is possibly too large!
  warnings.warn('The stepsize (%g) is possibly too large!' % h1[i])
model.params

OrderedDict([('shape', 0.050109518363545352),
             ('location', 3.8747498425529501),
             ('scale', 0.19804394476624812)])

model.plot_summary()
plt.show()
_images/User_guide_18_0.png