Skip to content

Plotting

pyextremes.plotting.extremes.plot_extremes(ts, extremes, extremes_method, extremes_type=None, block_size=None, threshold=None, r=None, figsize=(8, 5), ax=None)

Plot extreme events.

Parameters:

Name Type Description Default
ts Series

Time series from which extremes were extracted.

 required
extremes Series

Time series of extreme events.

 required
extremes_method str

Extreme value extraction method. Supported values: BM - Block Maxima POT - Peaks Over Threshold

 required
extremes_type str

Type of extremes, used only if extremes_method is 'POT' and threshold is not provided. high - extreme high values low - get low values

 None
block_size str or Timedelta

Block size, used only if extremes_method is 'BM'. If None (default), then calculated as median distance between extreme events.

 None
threshold float

Threshold, used only if extremes_method is 'POT'. If None (default), then is inferred from extremes as minimum if extremes_type is 'high' or maximum if extremes_type is 'low'.

 None
r pandas.Timedelta or value convertible to timedelta

Duration of window used to decluster the exceedances. See pandas.to_timedelta for more information. Used to show clusters. If None (default) then clusters are not shown. Clusters are shown only if both threshold and r were provided.

 None
figsize tuple

Figure size in inches in format (width, height). By default it is (8, 5).

 (8, 5)
ax Axes

Axes onto which extremes plot is drawn. If None (default), a new figure and axes objects are created.

 None

Returns:

Name Type Description
figure Figure

Figure object.
axes Axes

Axes object.
Source code in src/pyextremes/plotting/extremes.py 
 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
def plot_extremes(
    ts: pd.Series,
    extremes: pd.Series,
    extremes_method: Literal["BM", "POT"],
    extremes_type: Optional[Literal["high", "low"]] = None,
    block_size: Optional[Union[str, pd.Timedelta]] = None,
    threshold: Optional[float] = None,
    r: Optional[Union[pd.Timedelta, Any]] = None,
    figsize: Tuple[float, float] = (8, 5),
    ax: Optional[plt.Axes] = None,
) -> Tuple[plt.Figure, plt.Axes]:
    """
    Plot extreme events.

    Parameters
    ----------
    ts : pandas.Series
        Time series from which `extremes` were extracted.
    extremes : pandas.Series
        Time series of extreme events.
    extremes_method : str
        Extreme value extraction method.
        Supported values:
            BM - Block Maxima
            POT - Peaks Over Threshold
    extremes_type : str, optional
        Type of `extremes`, used only if `extremes_method` is 'POT'
        and `threshold` is not provided.
            high - extreme high values
            low - get low values
    block_size : str or pandas.Timedelta, optional
        Block size, used only if `extremes_method` is 'BM'.
        If None (default), then calculated as median distance between extreme events.
    threshold : float, optional
        Threshold, used only if `extremes_method` is 'POT'.
        If None (default), then is inferred from `extremes` as
        minimum if `extremes_type` is 'high' or maximum if `extremes_type` is 'low'.
    r : pandas.Timedelta or value convertible to timedelta, optional
        Duration of window used to decluster the exceedances.
        See pandas.to_timedelta for more information.
        Used to show clusters. If None (default) then clusters are not shown.
        Clusters are shown only if both `threshold` and `r` were provided.
    figsize : tuple, optional
        Figure size in inches in format (width, height).
        By default it is (8, 5).
    ax : matplotlib.axes._axes.Axes, optional
        Axes onto which extremes plot is drawn.
        If None (default), a new figure and axes objects are created.

    Returns
    -------
    figure : matplotlib.figure.Figure
        Figure object.
    axes : matplotlib.axes._axes.Axes
        Axes object.

    """
    if extremes_method not in ["BM", "POT"]:
        raise ValueError(
            f"invalid value in '{extremes_method}' for the 'extremes_method' argument"
        )

    if extremes_type not in ["high", "low"]:
        raise ValueError(
            f"invalid value in '{extremes_type}' for the 'extremes_type' argument"
        )

    with plt.rc_context(rc=pyextremes_rc):
        # Create figure
        if ax is None:
            fig, ax = plt.subplots(figsize=figsize, dpi=96)
        else:
            try:
                fig = ax.get_figure()
            except AttributeError as _error:
                raise TypeError(
                    f"invalid type in {type(ax)} for the 'ax' argument, "
                    f"must be matplotlib Axes object"
                ) from _error

        # Configure axes
        ax.grid(False)

        # Plot signal time series
        ax.plot(ts.index, ts.values, ls="-", color="#5199FF", lw=0.25, zorder=10)

        # Plot extreme events
        ax.scatter(
            extremes.index,
            extremes.values,
            s=20,
            lw=0.5,
            edgecolor="w",
            facecolor="#F85C50",
            zorder=20,
        )

        # Label the axes
        ax.set_xlabel(extremes.index.name or "date-time")
        ax.set_ylabel(extremes.name or "Extreme value")

        if extremes_method == "BM":
            # Infer 'block_size'
            if block_size is None:
                # Calculate 'block_size' as median of distances between extremes
                block_size = pd.to_timedelta(
                    np.quantile(
                        np.diff(extremes.index),
                        0.5,
                    )
                )
            else:
                if not isinstance(block_size, pd.Timedelta):
                    if isinstance(block_size, str):
                        block_size = pd.to_timedelta(block_size)
                    else:
                        raise TypeError(
                            f"invalid type in {type(block_size)} "
                            f"for the 'block_size' argument"
                        )

            # Plot block boundaries
            block_left_boundary = ts.index[0]
            while block_left_boundary < extremes.index.max() + block_size:
                ax.axvline(
                    block_left_boundary, ls="--", lw=0.5, color="#D1D3D4", zorder=5
                )
                block_left_boundary += block_size

        else:
            if threshold is None:
                if extremes_type is None:
                    raise TypeError(
                        "'extremes_type' argument must be provided "
                        "for 'extremes_method' being 'POT' "
                        "when 'threshold' is not provided"
                    )
                if extremes_type == "high":
                    threshold = extremes.min()
                else:
                    threshold = extremes.max()
            else:
                if r is not None:
                    # Plot clusters (only if both threshold and r are provided)
                    if extremes_type == "high":
                        exceedances = ts.loc[ts.values > threshold]
                    else:
                        exceedances = ts.loc[ts.values < threshold]
                    for cluster in _generate_clusters(exceedances=exceedances, r=r):
                        _plot_cluster(ax=ax, cluster=cluster)

            # Plot threshold line
            ax.axhline(threshold, ls="--", lw=1, color="#FF756B", zorder=15)

        fig.autofmt_xdate()

        return fig, ax

pyextremes.plotting.return_values.plot_return_values(observed_return_values, modeled_return_values, ax=None, figsize=(8, 5))

Plot return values and confidence intervals for given return periods.

Parameters:

Name Type Description Default
observed_return_values DataFrame

DataFrame with observed return values. First column must have extreme values. Must have 'return period' column.

 required
modeled_return_values DataFrame

DataFrame with modeled return values. Index has return periods. Must have the following columns: 'return value', 'lower ci', 'upper ci'.

 required
ax Axes

Axes onto which the return value plot is drawn. If None (default), a new figure and axes objects are created.

 None
figsize tuple

Figure size in inches in format (width, height). By default it is (8, 5).

 (8, 5)

Returns:

Name Type Description
figure Figure

Figure object.
axes Axes

Axes object.
Source code in src/pyextremes/plotting/return_values.py 
  9
 10
 11
 12
 13
 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
def plot_return_values(
    observed_return_values: pd.DataFrame,
    modeled_return_values: pd.DataFrame,
    ax: Optional[plt.Axes] = None,
    figsize: Tuple[float, float] = (8, 5),
) -> Tuple[plt.Figure, plt.Axes]:
    """
    Plot return values and confidence intervals for given return periods.

    Parameters
    ----------
    observed_return_values : pandas.DataFrame
        DataFrame with observed return values.
        First column must have extreme values.
        Must have 'return period' column.
    modeled_return_values : pandas.DataFrame
        DataFrame with modeled return values.
        Index has return periods.
        Must have the following columns: 'return value', 'lower ci', 'upper ci'.
    ax : matplotlib.axes._axes.Axes, optional
        Axes onto which the return value plot is drawn.
        If None (default), a new figure and axes objects are created.
    figsize : tuple, optional
        Figure size in inches in format (width, height).
        By default it is (8, 5).

    Returns
    -------
    figure : matplotlib.figure.Figure
        Figure object.
    axes : matplotlib.axes._axes.Axes
        Axes object.

    """
    # Validate the 'observed_return_values' argument
    if (
        len(observed_return_values.columns) < 2
        or observed_return_values.columns[0] == "return period"
        or "return period" not in observed_return_values.columns
    ):
        raise ValueError(
            f"'observed_return_values' argument "
            f"has invalid columns in {observed_return_values.columns}, "
            f"must have at least two columns and a 'return period' column "
            f"which is not the first column"
        )

    # Validate the 'modeled_return_values' argument
    if len(modeled_return_values.columns) < 3 or any(
        col not in modeled_return_values.columns
        for col in ["return value", "lower ci", "upper ci"]
    ):
        raise ValueError(
            f"'modeled_return_values' argument "
            f"has invalid columns in {modeled_return_values.columns}, "
            f"must have at least three columns and include columns: "
            f"'return value', 'lower ci', 'upper ci'"
        )

    with plt.rc_context(rc=pyextremes_rc):
        # Create figure
        if ax is None:
            fig, ax = plt.subplots(figsize=figsize, dpi=96)
        else:
            try:
                fig = ax.figure
            except AttributeError as _error:
                raise TypeError(
                    f"invalid type in {type(ax)} for the 'ax' argument, "
                    f"must be matplotlib Axes object"
                ) from _error

        # Configure axes
        ax.semilogx()
        ax.grid(True, which="both")
        ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x:,.0f}"))

        # Plot modeled confidence intervals
        for col in ["lower ci", "upper ci"]:
            ax.plot(
                modeled_return_values.index.values,
                modeled_return_values.loc[:, col].values,
                color="#5199FF",
                lw=1,
                ls="--",
                zorder=15,
            )
        ax.fill_between(
            modeled_return_values.index.values,
            modeled_return_values.loc[:, "lower ci"].values,
            modeled_return_values.loc[:, "upper ci"].values,
            facecolor="#5199FF",
            edgecolor="None",
            alpha=0.25,
            zorder=10,
        )

        # Plot observed extreme values
        ax.scatter(
            observed_return_values.loc[:, "return period"].values,
            observed_return_values.loc[:, observed_return_values.columns[0]].values,
            marker="o",
            s=20,
            lw=1,
            facecolor="k",
            edgecolor="w",
            zorder=20,
        )

        # Plot modeled return values
        ax.plot(
            modeled_return_values.index.values,
            modeled_return_values.loc[:, "return value"].values,
            color="#F85C50",
            lw=2,
            ls="-",
            zorder=25,
        )

        # Label axes
        ax.set_xlabel("Return period")
        ax.set_ylabel(observed_return_values.columns[0])

        return fig, ax

pyextremes.plotting.probability_plots.plot_probability(observed, theoretical, ax=None, figsize=(8, 8))

Plot a probability plot (QQ or PP).

Parameters:

Name Type Description Default
observed ndarray

Observed values.

 required
theoretical ndarray

Theoretical values.

 required
ax Axes

Axes onto which the probability plot is drawn. If None (default), a new figure and axes objects are created.

 None
figsize tuple

Figure size in inches in format (width, height). By default it is (8, 8).

 (8, 8)

Returns:

Name Type Description
figure Figure

Figure object.
axes Axes

Axes object.
Source code in src/pyextremes/plotting/probability_plots.py 
10
11
12
13
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
def plot_probability(
    observed: np.ndarray,
    theoretical: np.ndarray,
    ax: Optional[plt.Axes] = None,
    figsize: Tuple[float, float] = (8, 8),
) -> Tuple[plt.Figure, plt.Axes]:
    """
    Plot a probability plot (QQ or PP).

    Parameters
    ----------
    observed : numpy.ndarray
        Observed values.
    theoretical : numpy.ndarray
        Theoretical values.
    ax : matplotlib.axes._axes.Axes, optional
        Axes onto which the probability plot is drawn.
        If None (default), a new figure and axes objects are created.
    figsize : tuple, optional
        Figure size in inches in format (width, height).
        By default it is (8, 8).

    Returns
    -------
    figure : matplotlib.figure.Figure
        Figure object.
    axes : matplotlib.axes._axes.Axes
        Axes object.

    """
    with plt.rc_context(rc=pyextremes_rc):
        # Create figure
        if ax is None:
            fig, ax = plt.subplots(figsize=figsize, dpi=96)
        else:
            try:
                fig = ax.figure
            except AttributeError as _error:
                raise TypeError(
                    f"invalid type in {type(ax)} for the 'ax' argument, "
                    f"must be matplotlib Axes object"
                ) from _error

        # Configure axes
        ax.grid(False)

        # Plot scatter of observed and theoretical probabilities
        ax.scatter(
            theoretical,
            observed,
            marker="o",
            s=20,
            lw=0.75,
            facecolor="k",
            edgecolor="w",
            zorder=10,
        )

        # Plot a diagonal perfect-fit line
        min_value = min([min(ax.get_xlim()), min(ax.get_ylim())])
        max_value = max([max(ax.get_xlim()), max(ax.get_ylim())])
        ax.plot(
            [min_value, max_value],
            [min_value, max_value],
            color="#5199FF",
            lw=1,
            ls="--",
            zorder=5,
        )

        # Label axes
        ax.set_xlabel("Theoretical")
        ax.set_ylabel("Observed")

        # Calculate Pearson R statistic and show it in the figure
        pearsonr, p_value = scipy.stats.pearsonr(theoretical, observed)
        axes_range = max_value - min_value
        ax.text(
            x=min_value + 0.05 * axes_range,
            y=max_value - 0.05 * axes_range,
            s=f"$R^2={pearsonr:.3f}$\n$p={p_value:.3f}$",
            horizontalalignment="left",
            verticalalignment="top",
        )

        # Set axes limits
        ax.set_xlim(min_value, max_value)
        ax.set_ylim(min_value, max_value)

        return fig, ax

pyextremes.plotting.mcmc.plot_trace(trace, trace_map=None, burn_in=0, labels=None, figsize=None)

Plot a trace plot for a given MCMC sampler trace.

Parameters:

Name Type Description Default
trace ndarray

Array with MCMC sampler trace. Has a shape of (n_walkers, n_samples, n_parameters).

 required
trace_map tuple

Tuple with maximum aposteriori estimate of distribution parameters. If provided, MAP values are plotted as orange lines on top of the trace. If None (default) then MAP estimates are not plotted.

 None
burn_in int

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

 0
labels list of strings

Sequence of strings with parameter names, used to label axes. If None (default), then axes are labeled sequentially.

 None
figsize tuple

Figure size in inches. If None (default), then figure size is calculated automatically as 8 by 2 times number of parameters (trace.shape[2]).

 None

Returns:

Name Type Description
figure Figure

Figure object.
axes list

List with trace.shape[2] Axes objects.
Source code in src/pyextremes/plotting/mcmc.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
def plot_trace(
    trace: np.ndarray,
    trace_map: Optional[tuple] = None,
    burn_in: int = 0,
    labels: List[str] = None,
    figsize: Optional[Tuple[float, float]] = None,
) -> Tuple[plt.Figure, list]:
    """
    Plot a trace plot for a given MCMC sampler trace.

    Parameters
    ----------
    trace : numpy.ndarray
        Array with MCMC sampler trace.
        Has a shape of (n_walkers, n_samples, n_parameters).
    trace_map : tuple, optional
        Tuple with maximum aposteriori estimate of distribution parameters.
        If provided, MAP values are plotted as orange lines on top of the trace.
        If None (default) then MAP estimates are not plotted.
    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).
    labels : list of strings, optional
        Sequence of strings with parameter names, used to label axes.
        If None (default), then axes are labeled sequentially.
    figsize : tuple, optional
        Figure size in inches.
        If None (default), then figure size is calculated automatically
        as 8 by 2 times number of parameters (`trace.shape[2]`).

    Returns
    -------
    figure : matplotlib.figure.Figure
        Figure object.
    axes : list
        List with `trace.shape[2]` Axes objects.

    """
    # Parse the 'burn_in' argument
    if not isinstance(burn_in, int):
        raise TypeError(
            f"invalid type in {type(burn_in)} for the 'burn_in' argument, "
            f"must be integer"
        )
    if burn_in < 0:
        raise ValueError(
            f"invalid value in {burn_in} for the 'burn_in' argument, "
            f"must be a positive integer"
        )
    if burn_in >= trace.shape[1]:
        raise ValueError(
            f"'burn_in' value of {burn_in} exceeds number of samples {trace.shape[1]}"
        )

    # Calculate figure size
    n_parameters = trace.shape[2]
    figsize = figsize or (8, 2 * n_parameters)

    with plt.rc_context(rc=pyextremes_rc):
        # Create figure
        fig = plt.figure(figsize=figsize, dpi=96)

        # Create gridspec
        gs = matplotlib.gridspec.GridSpec(
            nrows=n_parameters,
            ncols=1,
            wspace=0.1,
            hspace=0.1,
            width_ratios=[1],
            height_ratios=np.full(shape=n_parameters, fill_value=1),
        )

        # Create and configure axes
        axes = [fig.add_subplot(gs[i, 0]) for i in range(n_parameters)]
        if labels is None:
            labels = [f"Parameter {i}" for i in range(n_parameters)]
        for i, ax in enumerate(axes):
            ax.grid(False)
            ax.set_ylabel(labels[i])
            if i == n_parameters - 1:
                ax.set_xlabel("Sample number")
                ax.xaxis.set_major_formatter(
                    plt.FuncFormatter(lambda x, _: f"{x:,.0f}"),
                )
            else:
                ax.tick_params(axis="x", which="both", labelbottom=False)

        # Plot the trace plots
        lines = np.full(
            shape=(trace.shape[0], trace.shape[1] - burn_in, 2),
            fill_value=np.nan,
            dtype=np.float64,
        )
        lines[:, :, 0] = np.arange(burn_in + 1, trace.shape[1] + 1, 1)
        for i, ax in enumerate(axes):
            lines[:, :, 1] = trace[:, burn_in:, i]
            line_collection = LineCollection(
                segments=lines.copy(),
                color="#231F20",
                lw=0.1,
                zorder=5,
            )
            ax.add_collection(line_collection)
            if trace_map is not None:
                ax.axhline(trace_map[i], color="#F85C50", lw=2, ls="--", zorder=10)
            ax.autoscale(enable=True, axis="both", tight=False)

        # Align y-labels
        fig.align_ylabels(axs=axes)

        return fig, axes

pyextremes.plotting.mcmc.plot_corner(trace, trace_map=None, burn_in=0, labels=None, levels=None, figsize=(8, 8))

Plot a corner plot for a given MCMC sampler trace.

Parameters:

Name Type Description Default
trace ndarray

Array with MCMC sampler trace. Has a shape of (n_walkers, n_samples, n_parameters).

 required
trace_map tuple

Tuple with maximum aposteriori estimate of distribution parameters. If provided, MAP values are plotted as orange lines. If None (default) then MAP estimates are not plotted.

 None
burn_in int

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

 0
labels array - like

Sequence of strings with parameter names, used to label axes. If None (default), then axes are labeled sequentially.

 None
levels int

Number of Gaussian KDE contours to plot. If None (default), then not shown.

 None
figsize tuple

Figure size in inches. By default it is (8, 8).

 (8, 8)

Returns:

Name Type Description
figure Figure

Figure object.
axes list

2D list with Axes objects of size N by N, where N is trace.shape[2]. Empty slots are represented by None. Axes are ordered from left to right top to bottom.
Source code in src/pyextremes/plotting/mcmc.py 
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
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
def plot_corner(
    trace: np.ndarray,
    trace_map: Optional[tuple] = None,
    burn_in: int = 0,
    labels: List[str] = None,
    levels: Optional[int] = None,
    figsize: Tuple[float, float] = (8, 8),
) -> Tuple[plt.Figure, list]:
    """
    Plot a corner plot for a given MCMC sampler trace.

    Parameters
    ----------
    trace : numpy.ndarray
        Array with MCMC sampler trace.
        Has a shape of (n_walkers, n_samples, n_parameters).
    trace_map : tuple, optional
        Tuple with maximum aposteriori estimate of distribution parameters.
        If provided, MAP values are plotted as orange lines.
        If None (default) then MAP estimates are not plotted.
    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).
    labels : array-like, optional
        Sequence of strings with parameter names, used to label axes.
        If None (default), then axes are labeled sequentially.
    levels : int, optional
        Number of Gaussian KDE contours to plot.
        If None (default), then not shown.
    figsize : tuple, optional
        Figure size in inches. By default it is (8, 8).

    Returns
    -------
    figure : matplotlib.figure.Figure
        Figure object.
    axes : list
        2D list with Axes objects of size N by N, where N is `trace.shape[2]`.
        Empty slots are represented by None. Axes are ordered from left to right
        top to bottom.

    """
    # Parse the 'burn_in' argument
    if not isinstance(burn_in, int):
        raise TypeError(
            f"invalid type in {type(burn_in)} for the 'burn_in' argument, "
            f"must be integer"
        )
    if burn_in < 0:
        raise ValueError(
            f"invalid value in {burn_in} for the 'burn_in' argument, "
            f"must be a positive integer"
        )
    if burn_in >= trace.shape[1]:
        raise ValueError(
            f"'burn_in' value of {burn_in} exceeds number of samples {trace.shape[1]}"
        )

    n_parameters = trace.shape[2]

    with plt.rc_context(rc=pyextremes_rc):
        # Create figure
        fig = plt.figure(figsize=figsize, dpi=96)

        # Create gridspec
        gs = matplotlib.gridspec.GridSpec(
            nrows=n_parameters,
            ncols=n_parameters,
            wspace=0.1,
            hspace=0.1,
            width_ratios=[1] * n_parameters,
            height_ratios=[1] * n_parameters,
        )

        # Create and configure axes
        axes = [[None] * n_parameters for _ in range(n_parameters)]
        for i in range(n_parameters):
            for j in range(n_parameters):
                # Create axes only for axes at or left of the main diagonal
                if i >= j:
                    ax = fig.add_subplot(gs[i, j])
                    ax.grid(False)
                    axes[i][j] = ax

                    if i == n_parameters - 1:
                        # Set x-axis labels for axes located in the first row
                        if labels is None:
                            ax.set_xlabel(f"Parameter {j}")
                        else:
                            ax.set_xlabel(labels[j])
                    else:
                        # Remove x-axis ticks for axes located above the bottom row
                        ax.tick_params(
                            axis="x", which="both", labelbottom=False, length=0
                        )

                    if j == 0:
                        # Set y-axis label for axes located in the first column
                        if labels is None:
                            ax.set_ylabel(f"Parameter {i}")
                        else:
                            ax.set_ylabel(labels[i])

                    if j != 0 or i == j == 0:
                        # Remove y-axis ticks for axes located right of the first column
                        # and for the first axes along the main diagonal
                        ax.tick_params(
                            axis="y", which="both", labelleft=False, length=0
                        )

                    if i == j:
                        # Plot histogram
                        parameter_samples = trace[:, burn_in:, i].flatten()
                        ax.hist(
                            parameter_samples,
                            bins=np.histogram_bin_edges(
                                a=parameter_samples, bins="auto"
                            ),
                            density=True,
                            histtype="step",
                            edgecolor="#231F20",
                            lw=0.5,
                            zorder=5,
                        )
                        if trace_map is not None:
                            ax.axvline(
                                trace_map[i], color="#F85C50", lw=1, ls="--", zorder=10
                            )
                    else:
                        # Plot scatter plot
                        parameter_i = trace[:, burn_in:, i].flatten()
                        parameter_j = trace[:, burn_in:, j].flatten()
                        ax.scatter(
                            parameter_j,
                            parameter_i,
                            marker="o",
                            s=2,
                            alpha=0.1,
                            facecolor="#231F20",
                            edgecolor="None",
                            lw=0,
                            zorder=5,
                        )

                        # Plot trace map lines
                        if trace_map is not None:
                            ax.axvline(
                                trace_map[j], color="#F85C50", lw=1, ls="--", zorder=15
                            )
                            ax.axhline(
                                trace_map[i], color="#F85C50", lw=1, ls="--", zorder=15
                            )

                        # Plot Gaussian KDE contour
                        if levels is not None:
                            kernel = scipy.stats.gaussian_kde(
                                np.vstack([parameter_j, parameter_i])
                            )
                            xx, yy = np.meshgrid(
                                np.linspace(parameter_j.min(), parameter_j.max(), 100),
                                np.linspace(parameter_i.min(), parameter_i.max(), 100),
                            )
                            zz = np.reshape(
                                kernel(np.vstack([xx.ravel(), yy.ravel()])).transpose(),
                                xx.shape,
                            )
                            ax.contour(
                                xx,
                                yy,
                                zz,
                                colors="w",
                                levels=levels,
                                linewidths=1,
                                linestyles="-",
                                zorder=10,
                            )

                    # Set axes limits
                    parameter_i = trace[:, burn_in:, i].flatten()
                    parameter_j = trace[:, burn_in:, j].flatten()
                    xlim = (parameter_j.min(), parameter_j.max())
                    ylim = (parameter_i.min(), parameter_i.max())
                    ax.set_xlim(xlim)
                    if i != j:
                        ax.set_ylim(ylim)

        # Align labels
        fig.align_labels()

        return fig, axes