Skip to content

Models

pyextremes.models.models.get_model(model, extremes, distribution, distribution_kwargs=None, **kwargs) 

get_model(
    model: Literal["MLE"],
    extremes: pd.Series,
    distribution: Union[str, scipy.stats.rv_continuous],
    distribution_kwargs: Optional[dict] = None,
) -> MLE

get_model(
    model: Literal["MOM"],
    extremes: pd.Series,
    distribution: Union[str, scipy.stats.rv_continuous],
    distribution_kwargs: Optional[dict] = None,
) -> MOM

get_model(
    model: Literal["Lmoments"],
    extremes: pd.Series,
    distribution: Union[str, scipy.stats.rv_continuous],
    distribution_kwargs: Optional[dict] = None,
) -> Lmoments

get_model(
    model: Literal["Emcee"],
    extremes: pd.Series,
    distribution: Union[str, scipy.stats.rv_continuous],
    distribution_kwargs: Optional[dict] = None,
    *,
    n_walkers: int = 100,
    n_samples: int = 500,
    progress: bool = False,
) -> Emcee

Get distribution fitting model and fit it to given extreme values.

Parameters:

Name Type Description Default
model str

Name of model. Supported models: MLE - Maximum Likelihood Estimate (MLE) model. Based on 'scipy' package (scipy.stats.rv_continuous.fit). Emcee - Markov Chain Monte Carlo (MCMC) model. Based on 'emcee' package by Daniel Foreman-Mackey. Lmoments - L-moments (Lmoments) model. Based on 'lmoments3' package, https://github.com/Ouranosinc/lmoments3 MOM - Method of Moments (MOM) model. Based on 'scipy' package (scipy.stats.rv_continuous.fit).

 required
extremes Series

Time series of extreme events.

 required
distribution str or rv_continuous

Distribution name compatible with scipy.stats or a subclass of scipy.stats.rv_continuous. See https://docs.scipy.org/doc/scipy/reference/stats.html

 required
distribution_kwargs dict

Special keyword arguments, passed to the .fit method of the distribution. These keyword arguments represent parameters to be held fixed. Names of parameters to be fixed must have 'f' prefixes. Valid parameters: - shape(s): 'fc', e.g. fc=0 - location: 'floc', e.g. floc=0 - scale: 'fscale', e.g. fscale=1 By default, no parameters are fixed. See documentation of a specific scipy.stats distribution for names of available parameters.

 None
kwargs

Keyword arguments passed to a model .fit method. MLE model: MLE model takes no additional arguments. Emcee model: n_walkers : int, optional The number of walkers in the ensemble (default=100). n_samples : int, optional The number of steps to run (default=500). progress : bool or str, optional If True, a progress bar will be shown as the sampler progresses. If a string, will select a specific tqdm progress bar. Most notable is 'notebook', which shows a progress bar suitable for Jupyter notebooks. If False (default), no progress bar will be shown. This progress bar is a part of the emcee package.

 {}

Returns:

Name Type Description
model MLE, Emcee, Lmoments, or MOM
Distribution fitting model fitted to the `extremes`.
Source code in src/pyextremes/models/models.py 
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
def get_model(
    model: Literal["MLE", "Emcee", "Lmoments", "MOM"],
    extremes: pd.Series,
    distribution: Union[str, scipy.stats.rv_continuous],
    distribution_kwargs: Optional[dict] = None,
    **kwargs,
) -> Union[MLE, Emcee, Lmoments, MOM]:
    """
    Get distribution fitting model and fit it to given extreme values.

    Parameters
    ----------
    model : str
        Name of model.
        Supported models:
            MLE - Maximum Likelihood Estimate (MLE) model.
                Based on 'scipy' package (scipy.stats.rv_continuous.fit).
            Emcee - Markov Chain Monte Carlo (MCMC) model.
                Based on 'emcee' package by Daniel Foreman-Mackey.
            Lmoments - L-moments (Lmoments) model.
                Based on 'lmoments3' package, https://github.com/Ouranosinc/lmoments3
            MOM - Method of Moments (MOM) model.
                Based on 'scipy' package (scipy.stats.rv_continuous.fit).
    extremes : pandas.Series
        Time series of extreme events.
    distribution : str or scipy.stats.rv_continuous
        Distribution name compatible with scipy.stats
        or a subclass of scipy.stats.rv_continuous.
        See https://docs.scipy.org/doc/scipy/reference/stats.html
    distribution_kwargs : dict, optional
        Special keyword arguments, passed to the `.fit` method of the distribution.
        These keyword arguments represent parameters to be held fixed.
        Names of parameters to be fixed must have 'f' prefixes. Valid parameters:
            - shape(s): 'fc', e.g. fc=0
            - location: 'floc', e.g. floc=0
            - scale: 'fscale', e.g. fscale=1
        By default, no parameters are fixed.
        See documentation of a specific scipy.stats distribution
        for names of available parameters.
    kwargs
        Keyword arguments passed to a model .fit method.
        MLE model:
            MLE model takes no additional arguments.
        Emcee model:
            n_walkers : int, optional
                The number of walkers in the ensemble (default=100).
            n_samples : int, optional
                The number of steps to run (default=500).
            progress : bool or str, optional
                If True, a progress bar will be shown as the sampler progresses.
                If a string, will select a specific tqdm progress bar.
                Most notable is 'notebook', which shows a progress bar
                suitable for Jupyter notebooks.
                If False (default), no progress bar will be shown.
                This progress bar is a part of the `emcee` package.

    Returns
    -------
    model : MLE, Emcee, Lmoments, or MOM
    Distribution fitting model fitted to the `extremes`.

    """
    distribution_model_kwargs = {
        "extremes": extremes,
        "distribution": distribution,
        "distribution_kwargs": distribution_kwargs,
        **kwargs,
    }

    if model == "MLE":
        return MLE(**distribution_model_kwargs)
    if model == "Emcee":
        return Emcee(**distribution_model_kwargs)
    if model == "MOM":
        return MOM(**distribution_model_kwargs)
    if model == "Lmoments":
        return Lmoments(**distribution_model_kwargs)
    raise ValueError(
        f"invalid value in '{model}' for the 'model' argument, "
        f"available models: 'MLE', 'Emcee', 'Lmoments', 'MOM'"
    )

pyextremes.models.model_base.AbstractModelBaseClass

Bases: ABC

Source code in src/pyextremes/models/model_base.py 
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
class AbstractModelBaseClass(abc.ABC):
    def __init__(
        self,
        extremes: pd.Series,
        distribution: typing.Union[str, scipy.stats.rv_continuous],
        distribution_kwargs: typing.Optional[dict] = None,
        **kwargs,
    ) -> None:
        """
        Initialize the model.

        Parameters
        ----------
        extremes : pandas.Series
            Time series of extreme events.
        distribution : str or scipy.stats.rv_continuous
            Distribution name compatible with scipy.stats
            or a subclass of scipy.stats.rv_continuous.
            See https://docs.scipy.org/doc/scipy/reference/stats.html
        distribution_kwargs : dict, optional
            Special keyword arguments, passed to the `.fit` method of the distribution.
            These keyword arguments represent parameters to be held fixed.
            Names of parameters to be fixed must have 'f' prefixes. Valid parameters:
                - shape(s): 'fc', e.g. fc=0
                - location: 'floc', e.g. floc=0
                - scale: 'fscale', e.g. fscale=1
            By default, no parameters are fixed.
            See documentation of a specific scipy.stats distribution
            for names of available parameters.
            Lmoments does not allow fixed parameters.
        kwargs
            Keyword arguments passed to a model .fit method.
            MLE model:
                MLE model takes no additional arguments.
            MOM model:
                MOM model takes no additional arguments.
            Lmoments model:
                Lmoments model takes no additional arguments.
            Emcee model:
                n_walkers : int, optional
                    The number of walkers in the ensemble (default=100).
                n_samples : int, optional
                    The number of steps to run (default=500).
                progress : bool or str, optional
                    If True, a progress bar will be shown as the sampler progresses.
                    If a string, will select a specific tqdm progress bar.
                    Most notable is 'notebook', which shows a progress bar
                    suitable for Jupyter notebooks.
                    If False (default), no progress bar will be shown.
                    This progress bar is a part of the `emcee` package.

        """
        self.extremes = extremes

        # Declare extreme value distribution
        distribution_kwargs = distribution_kwargs or {}
        self.distribution = Distribution(
            extremes=self.extremes,
            distribution=distribution,
            **distribution_kwargs,
        )

        # Fit the distribution to extremes
        self._fit_parameters: typing.Optional[dict] = None
        self._trace: typing.Optional[np.ndarray] = None
        self.fit(**kwargs)

        # Initialize 'return_value_cache'
        self.return_value_cache: typing.Dict[tuple, tuple] = {}

    @property
    @abc.abstractmethod
    def name(self) -> str:
        """Return model name."""
        raise NotImplementedError

    @abc.abstractmethod
    def fit(self, **kwargs) -> None:
        """
        Set values for self.fit_parameters and self.trace.

        self.trace is set only for MCMC-like models.
        self.fit_parameters is a dictionary with {parameter_name: value},
        e.g. {'c': 0.1, 'loc': -7, 'scale': 0.3}
        self.trace is a numpy.ndarray with shape of
        (n_walkers, n_samples, n_free_parameters)

        """
        raise NotImplementedError

    @property
    def fit_parameters(self) -> typing.Dict[str, float]:
        if self._fit_parameters is None:
            raise AssertionError
        else:
            return self._fit_parameters

    @property
    def trace(self) -> np.ndarray:
        if self._trace is None:
            raise TypeError(f"trace property is not applicable for '{self.name}' model")
        else:
            return self._trace

    @property
    def loglikelihood(self) -> float:
        return np.sum(self.logpdf(x=self.extremes.values))

    @property
    def AIC(self) -> float:  # noqa: N802
        """
        Return corrected Akaike Information Criterion (AIC) of the model.

        Smaller AIC value corresponds to better model.
        This formula scales well for small sample sizes.
        See https://en.wikipedia.org/wiki/Akaike_information_criterion

        """
        k = self.distribution.number_of_parameters
        n = len(self.extremes)
        aic = 2 * k - 2 * self.loglikelihood
        correction = (2 * k**2 + 2 * k) / (n - k - 1)
        return aic + correction

    @abc.abstractmethod
    def get_return_value(
        self,
        exceedance_probability: float,
        alpha: typing.Optional[float] = None,
        **kwargs,
    ) -> tuple:
        """
        Calculate return value and confidence interval bounds.

        Parameters
        ----------
        exceedance_probability : array-like
            Exceedance probability or 1D array of exceedance probabilities.
            Each exceedance probability must be in the [0, 1) range.
        alpha : float, optional
            Width of confidence interval (0, 1).
            If None (default), return None
            for upper and lower confidence interval bounds.
        kwargs
            Model-specific keyword arguments.
            If alpha is None, keyword arguments are ignored
            (error still raised for unrecognized arguments).
            MLE model:
                n_samples : int, optional
                    Number of bootstrap samples used to estimate
                    confidence interval bounds (default=100).
            Emcee model:
                burn_in : int
                    Burn-in value (number of first steps to discard for each walker).

        Returns
        -------
        return_value : array-like
            Return values.
        ci_lower : array-like
            Lower confidence interval bounds.
        ci_upper : array-like
            Upper confidence interval bounds.

        """
        raise NotImplementedError

    def _get_prop(self, prop: str, x):
        return self.distribution.get_prop(
            prop=prop, x=x, free_parameters=self.fit_parameters
        )

    def pdf(self, x):
        return self._get_prop(prop="pdf", x=x)

    def logpdf(self, x):
        return self._get_prop(prop="logpdf", x=x)

    def cdf(self, x):
        return self._get_prop(prop="cdf", x=x)

    def ppf(self, x):
        return self._get_prop(prop="ppf", x=x)

    def isf(self, x):
        return self._get_prop(prop="isf", x=x)

AIC property

Return corrected Akaike Information Criterion (AIC) of the model.

Smaller AIC value corresponds to better model. This formula scales well for small sample sizes. See https://en.wikipedia.org/wiki/Akaike_information_criterion

name abstractmethod property

Return model name.

__init__(extremes, distribution, distribution_kwargs=None, **kwargs)

Initialize the model.

Parameters:

Name Type Description Default
extremes Series

Time series of extreme events.

 required
distribution str or rv_continuous

Distribution name compatible with scipy.stats or a subclass of scipy.stats.rv_continuous. See https://docs.scipy.org/doc/scipy/reference/stats.html

 required
distribution_kwargs dict

Special keyword arguments, passed to the .fit method of the distribution. These keyword arguments represent parameters to be held fixed. Names of parameters to be fixed must have 'f' prefixes. Valid parameters: - shape(s): 'fc', e.g. fc=0 - location: 'floc', e.g. floc=0 - scale: 'fscale', e.g. fscale=1 By default, no parameters are fixed. See documentation of a specific scipy.stats distribution for names of available parameters. Lmoments does not allow fixed parameters.

 None
kwargs

Keyword arguments passed to a model .fit method. MLE model: MLE model takes no additional arguments. MOM model: MOM model takes no additional arguments. Lmoments model: Lmoments model takes no additional arguments. Emcee model: n_walkers : int, optional The number of walkers in the ensemble (default=100). n_samples : int, optional The number of steps to run (default=500). progress : bool or str, optional If True, a progress bar will be shown as the sampler progresses. If a string, will select a specific tqdm progress bar. Most notable is 'notebook', which shows a progress bar suitable for Jupyter notebooks. If False (default), no progress bar will be shown. This progress bar is a part of the emcee package.

 {}
Source code in src/pyextremes/models/model_base.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
def __init__(
    self,
    extremes: pd.Series,
    distribution: typing.Union[str, scipy.stats.rv_continuous],
    distribution_kwargs: typing.Optional[dict] = None,
    **kwargs,
) -> None:
    """
    Initialize the model.

    Parameters
    ----------
    extremes : pandas.Series
        Time series of extreme events.
    distribution : str or scipy.stats.rv_continuous
        Distribution name compatible with scipy.stats
        or a subclass of scipy.stats.rv_continuous.
        See https://docs.scipy.org/doc/scipy/reference/stats.html
    distribution_kwargs : dict, optional
        Special keyword arguments, passed to the `.fit` method of the distribution.
        These keyword arguments represent parameters to be held fixed.
        Names of parameters to be fixed must have 'f' prefixes. Valid parameters:
            - shape(s): 'fc', e.g. fc=0
            - location: 'floc', e.g. floc=0
            - scale: 'fscale', e.g. fscale=1
        By default, no parameters are fixed.
        See documentation of a specific scipy.stats distribution
        for names of available parameters.
        Lmoments does not allow fixed parameters.
    kwargs
        Keyword arguments passed to a model .fit method.
        MLE model:
            MLE model takes no additional arguments.
        MOM model:
            MOM model takes no additional arguments.
        Lmoments model:
            Lmoments model takes no additional arguments.
        Emcee model:
            n_walkers : int, optional
                The number of walkers in the ensemble (default=100).
            n_samples : int, optional
                The number of steps to run (default=500).
            progress : bool or str, optional
                If True, a progress bar will be shown as the sampler progresses.
                If a string, will select a specific tqdm progress bar.
                Most notable is 'notebook', which shows a progress bar
                suitable for Jupyter notebooks.
                If False (default), no progress bar will be shown.
                This progress bar is a part of the `emcee` package.

    """
    self.extremes = extremes

    # Declare extreme value distribution
    distribution_kwargs = distribution_kwargs or {}
    self.distribution = Distribution(
        extremes=self.extremes,
        distribution=distribution,
        **distribution_kwargs,
    )

    # Fit the distribution to extremes
    self._fit_parameters: typing.Optional[dict] = None
    self._trace: typing.Optional[np.ndarray] = None
    self.fit(**kwargs)

    # Initialize 'return_value_cache'
    self.return_value_cache: typing.Dict[tuple, tuple] = {}

fit(**kwargs) abstractmethod

Set values for self.fit_parameters and self.trace.

self.trace is set only for MCMC-like models. self.fit_parameters is a dictionary with {parameter_name: value}, e.g. {'c': 0.1, 'loc': -7, 'scale': 0.3} self.trace is a numpy.ndarray with shape of (n_walkers, n_samples, n_free_parameters)

Source code in src/pyextremes/models/model_base.py 
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
@abc.abstractmethod
def fit(self, **kwargs) -> None:
    """
    Set values for self.fit_parameters and self.trace.

    self.trace is set only for MCMC-like models.
    self.fit_parameters is a dictionary with {parameter_name: value},
    e.g. {'c': 0.1, 'loc': -7, 'scale': 0.3}
    self.trace is a numpy.ndarray with shape of
    (n_walkers, n_samples, n_free_parameters)

    """
    raise NotImplementedError

get_return_value(exceedance_probability, alpha=None, **kwargs) abstractmethod

Calculate return value and confidence interval bounds.

Parameters:

Name Type Description Default
exceedance_probability array - like

Exceedance probability or 1D array of exceedance probabilities. Each exceedance probability must be in the [0, 1) range.

 required
alpha float

Width of confidence interval (0, 1). If None (default), return None for upper and lower confidence interval bounds.

 None
kwargs

Model-specific keyword arguments. If alpha is None, keyword arguments are ignored (error still raised for unrecognized arguments). MLE model: n_samples : int, optional Number of bootstrap samples used to estimate confidence interval bounds (default=100). Emcee model: burn_in : int Burn-in value (number of first steps to discard for each walker).

 {}

Returns:

Name Type Description
return_value array - like

Return values.
ci_lower array - like

Lower confidence interval bounds.
ci_upper array - like

Upper confidence interval bounds.
Source code in src/pyextremes/models/model_base.py 
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
@abc.abstractmethod
def get_return_value(
    self,
    exceedance_probability: float,
    alpha: typing.Optional[float] = None,
    **kwargs,
) -> tuple:
    """
    Calculate return value and confidence interval bounds.

    Parameters
    ----------
    exceedance_probability : array-like
        Exceedance probability or 1D array of exceedance probabilities.
        Each exceedance probability must be in the [0, 1) range.
    alpha : float, optional
        Width of confidence interval (0, 1).
        If None (default), return None
        for upper and lower confidence interval bounds.
    kwargs
        Model-specific keyword arguments.
        If alpha is None, keyword arguments are ignored
        (error still raised for unrecognized arguments).
        MLE model:
            n_samples : int, optional
                Number of bootstrap samples used to estimate
                confidence interval bounds (default=100).
        Emcee model:
            burn_in : int
                Burn-in value (number of first steps to discard for each walker).

    Returns
    -------
    return_value : array-like
        Return values.
    ci_lower : array-like
        Lower confidence interval bounds.
    ci_upper : array-like
        Upper confidence interval bounds.

    """
    raise NotImplementedError

pyextremes.models.model_emcee.Emcee

Bases: AbstractModelBaseClass

Markov Chain Monte Carlo (MCMC) model.

Built around the 'emcee' package by Daniel Foreman-Mackey

Source code in src/pyextremes/models/model_emcee.py 
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
class Emcee(AbstractModelBaseClass):
    """
    Markov Chain Monte Carlo (MCMC) model.

    Built around the 'emcee' package by Daniel Foreman-Mackey

    """

    def __init__(
        self,
        extremes: pd.Series,
        distribution: typing.Union[str, scipy.stats.rv_continuous],
        distribution_kwargs: typing.Optional[dict] = None,
        n_walkers: int = 100,
        n_samples: int = 500,
        progress: bool = False,
    ) -> None:
        super().__init__(
            extremes=extremes,
            distribution=distribution,
            distribution_kwargs=distribution_kwargs,
            n_walkers=n_walkers,
            n_samples=n_samples,
            progress=progress,
        )
        self.n_walkers = n_walkers
        self.n_samples = n_samples

    @property
    def name(self) -> str:
        return "Emcee"

    def fit(self, **kwargs) -> None:
        # Parse kwargs
        n_walkers: int = kwargs.pop("n_walkers")
        n_samples: int = kwargs.pop("n_samples")
        progress: bool = kwargs.pop("progress")
        if len(kwargs) != 0:
            raise TypeError(
                f"unrecognized arguments passed in: {', '.join(kwargs.keys())}"
            )

        # Declare Emcee ensemble sampler
        sampler = emcee.EnsembleSampler(
            nwalkers=n_walkers,
            ndim=self.distribution.number_of_parameters,
            log_prob_fn=self.distribution.log_probability,
        )

        # Run the ensemble sampler
        logger.debug(
            "running ensemble sampler with %d walkers and %d samples",
            n_walkers,
            n_samples,
        )
        with warnings.catch_warnings():
            warnings.simplefilter(action="ignore", category=RuntimeWarning)
            sampler.run_mcmc(
                initial_state=self.distribution.get_initial_state(n_walkers=n_walkers),
                nsteps=n_samples,
                progress=progress,
            )
        logger.debug(
            "finished run for ensemble sampler with %d walkers and %d samples",
            n_walkers,
            n_samples,
        )

        # Extract ensemble sampler chain
        self._trace: np.ndarray = sampler.get_chain().transpose((1, 0, 2))

        # Calculate fit parameters as MAP of distribution parameters
        kernel = scipy.stats.gaussian_kde(np.vstack(self._trace).transpose())

        def kde_func(x):
            return -kernel(x)[0]

        fit_parameters = self._trace.mean(axis=(0, 1))
        solution = scipy.optimize.minimize(
            kde_func,
            x0=fit_parameters,
            method="Nelder-Mead",
        )
        if solution.success:
            fit_parameters = solution.x
        else:  # pragma: no cover
            warnings.warn(
                message=(
                    "cannot calculate MAP using Gaussian KDE, "
                    "setting fit parameters as mean"
                ),
                category=RuntimeWarning,
            )
        self._fit_parameters = dict(
            zip(self.distribution.free_parameters, fit_parameters)
        )

    @property
    def trace_map(self) -> tuple:
        return tuple(
            self.fit_parameters[parameter]
            for parameter in self.distribution.free_parameters
        )

    def get_return_value(
        self, exceedance_probability, alpha: typing.Optional[float] = None, **kwargs
    ) -> tuple:
        """
        Calculate return value and confidence interval bounds.

        Parameters
        ----------
        exceedance_probability : array-like
            Exceedance probability or 1D array of exceedance probabilities.
            Each exceedance probability must be in the [0, 1) range.
        alpha : float, optional
            Width of confidence interval (0, 1).
            If None (default), return None
            for upper and lower confidence interval bounds.
        kwargs
            burn_in : int, optional
                Burn-in value (number of first steps to discard for each walker).
                By default it is 0 (no values are discarded).

        Returns
        -------
        return_value : array-like
            Return values.
        ci_lower : array-like
            Lower confidence interval bounds.
        ci_upper : array-like
            Upper confidence interval bounds.

        """
        # Parse 'kwargs'
        burn_in = kwargs.pop("burn_in", 0)
        if len(kwargs) != 0:
            raise TypeError(
                f"unrecognized arguments passed in: {', '.join(kwargs.keys())}"
            )

        # Convert 'exceedance_probability' to ndarray
        exceedance_probability = np.asarray(
            a=exceedance_probability, dtype=np.float64
        ).copy()
        if exceedance_probability.ndim == 0:
            exceedance_probability = exceedance_probability[np.newaxis]
        if exceedance_probability.ndim != 1:
            raise ValueError(
                f"invalid shape in {exceedance_probability.shape} "
                f"for the 'exceedance_probability' argument, must be 1D array"
            )

        # Calculate return values
        return_value = np.full(
            shape=exceedance_probability.shape, fill_value=np.nan, dtype=np.float64
        )
        ci_lower = return_value.copy()
        ci_upper = return_value.copy()
        for i, ep in enumerate(exceedance_probability):
            key: typing.Tuple[float, typing.Optional[float], int] = (
                ep,
                alpha,
                burn_in,
            )
            try:
                # Try to fetch pre-calculated values from cache
                rv, cil, ciu = self.return_value_cache[key]
                logger.debug(
                    "fetched return value for %s from cache as (%s, %s, %s)",
                    key,
                    rv,
                    cil,
                    ciu,
                )
            except KeyError:
                # Value not in cache - calculate new return value
                rv = self.distribution.distribution.isf(
                    q=ep,
                    **self.fit_parameters,
                    **self.distribution._fixed_parameters,
                )

                # Calculate confidence intervals
                if alpha is None:
                    cil = None
                    ciu = None
                else:
                    # Calculate confidence intervals
                    rv_sample = self.distribution.get_prop(
                        prop="isf",
                        x=ep,
                        free_parameters=np.vstack(self.trace[:, burn_in:, :]),
                    )
                    cil, ciu = np.quantile(
                        a=rv_sample, q=[(1 - alpha) / 2, (1 + alpha) / 2]
                    )

                # Add calculated return value and intervals to cache
                self.return_value_cache[key] = (rv, cil, ciu)
                logger.debug(
                    "calculated return value for %s as (%s, %s, %s)",
                    key,
                    rv,
                    cil,
                    ciu,
                )

            return_value[i] = rv
            ci_lower[i] = cil
            ci_upper[i] = ciu

        # Return results
        if len(return_value) == 1:
            return return_value[0], ci_lower[0], ci_upper[0]
        else:
            return return_value, ci_lower, ci_upper

    def __repr__(self) -> str:
        free_parameters = ", ".join(
            [
                f"{parameter}={self.fit_parameters[parameter]:.3f}"
                for parameter in self.distribution.free_parameters
            ]
        )

        fixed_parameters = ", ".join(
            [
                f"{key}={value:.3f}"
                for key, value in self.distribution.fixed_parameters.items()
            ]
        )
        if fixed_parameters == "":
            fixed_parameters = "all parameters are free"

        summary = [
            "Emcee model",
            "",
            f"free parameters: {free_parameters}",
            f"fixed parameters: {fixed_parameters}",
            f"AIC: {self.AIC:.3f}",
            f"loglikelihood: {self.loglikelihood:.3f}",
            f"number of walkers: {self.n_walkers:d}",
            f"number of samples: {self.n_samples:d}",
            f"return value cache size: {len(self.return_value_cache):,d}",
        ]

        longest_row = max(map(len, summary))
        summary[1] = "-" * longest_row
        summary.append(summary[1])
        summary[0] = " " * ((longest_row - len(summary[0])) // 2) + summary[0]
        for i, row in enumerate(summary):
            summary[i] += " " * (longest_row - len(row))

        return "\n".join(summary)

get_return_value(exceedance_probability, alpha=None, **kwargs)

Calculate return value and confidence interval bounds.

Parameters:

Name Type Description Default
exceedance_probability array - like

Exceedance probability or 1D array of exceedance probabilities. Each exceedance probability must be in the [0, 1) range.

 required
alpha float

Width of confidence interval (0, 1). If None (default), return None for upper and lower confidence interval bounds.

 None
kwargs

burn_in : int, optional Burn-in value (number of first steps to discard for each walker). By default it is 0 (no values are discarded).

 {}

Returns:

Name Type Description
return_value array - like

Return values.
ci_lower array - like

Lower confidence interval bounds.
ci_upper array - like

Upper confidence interval bounds.
Source code in src/pyextremes/models/model_emcee.py 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
def get_return_value(
    self, exceedance_probability, alpha: typing.Optional[float] = None, **kwargs
) -> tuple:
    """
    Calculate return value and confidence interval bounds.

    Parameters
    ----------
    exceedance_probability : array-like
        Exceedance probability or 1D array of exceedance probabilities.
        Each exceedance probability must be in the [0, 1) range.
    alpha : float, optional
        Width of confidence interval (0, 1).
        If None (default), return None
        for upper and lower confidence interval bounds.
    kwargs
        burn_in : int, optional
            Burn-in value (number of first steps to discard for each walker).
            By default it is 0 (no values are discarded).

    Returns
    -------
    return_value : array-like
        Return values.
    ci_lower : array-like
        Lower confidence interval bounds.
    ci_upper : array-like
        Upper confidence interval bounds.

    """
    # Parse 'kwargs'
    burn_in = kwargs.pop("burn_in", 0)
    if len(kwargs) != 0:
        raise TypeError(
            f"unrecognized arguments passed in: {', '.join(kwargs.keys())}"
        )

    # Convert 'exceedance_probability' to ndarray
    exceedance_probability = np.asarray(
        a=exceedance_probability, dtype=np.float64
    ).copy()
    if exceedance_probability.ndim == 0:
        exceedance_probability = exceedance_probability[np.newaxis]
    if exceedance_probability.ndim != 1:
        raise ValueError(
            f"invalid shape in {exceedance_probability.shape} "
            f"for the 'exceedance_probability' argument, must be 1D array"
        )

    # Calculate return values
    return_value = np.full(
        shape=exceedance_probability.shape, fill_value=np.nan, dtype=np.float64
    )
    ci_lower = return_value.copy()
    ci_upper = return_value.copy()
    for i, ep in enumerate(exceedance_probability):
        key: typing.Tuple[float, typing.Optional[float], int] = (
            ep,
            alpha,
            burn_in,
        )
        try:
            # Try to fetch pre-calculated values from cache
            rv, cil, ciu = self.return_value_cache[key]
            logger.debug(
                "fetched return value for %s from cache as (%s, %s, %s)",
                key,
                rv,
                cil,
                ciu,
            )
        except KeyError:
            # Value not in cache - calculate new return value
            rv = self.distribution.distribution.isf(
                q=ep,
                **self.fit_parameters,
                **self.distribution._fixed_parameters,
            )

            # Calculate confidence intervals
            if alpha is None:
                cil = None
                ciu = None
            else:
                # Calculate confidence intervals
                rv_sample = self.distribution.get_prop(
                    prop="isf",
                    x=ep,
                    free_parameters=np.vstack(self.trace[:, burn_in:, :]),
                )
                cil, ciu = np.quantile(
                    a=rv_sample, q=[(1 - alpha) / 2, (1 + alpha) / 2]
                )

            # Add calculated return value and intervals to cache
            self.return_value_cache[key] = (rv, cil, ciu)
            logger.debug(
                "calculated return value for %s as (%s, %s, %s)",
                key,
                rv,
                cil,
                ciu,
            )

        return_value[i] = rv
        ci_lower[i] = cil
        ci_upper[i] = ciu

    # Return results
    if len(return_value) == 1:
        return return_value[0], ci_lower[0], ci_upper[0]
    else:
        return return_value, ci_lower, ci_upper

pyextremes.models.model_mle.MLE

Bases: ScipyModel

Maximum Likelihood Estimate (MLE) model.

Built around the scipy.stats.rv_continuous.fit method.

Source code in src/pyextremes/models/model_mle.py 
 8
 9
10
11
12
13
14
15
16
17
18
class MLE(ScipyModel):
    """
    Maximum Likelihood Estimate (MLE) model.

    Built around the scipy.stats.rv_continuous.fit method.

    """

    @property
    def name(self) -> str:
        return "MLE"