# GEV responses
In this tutorial, we illustrate how to set up a distributional
regression model with the generalized extreme value distribution as a
response distribution. First, we simulate some data in R:
- The location parameter ($\mu$) is a function of an intercept and a
non-linear covariate effect.
- The scale parameter ($\sigma$) is a function of an intercept and a
linear effect and uses a log-link.
- The shape or concentration parameter ($\xi$) is a function of an
intercept and a linear effect.
After simulating the data, we can configure the model with a single call
to the `rliesel::liesel()` function.
``` r
library(rliesel)
```
Please make sure you are using a virtual or conda environment with Liesel installed, e.g. using `reticulate::use_virtualenv()` or `reticulate::use_condaenv()`. See `vignette("versions", "reticulate")`.
After setting the environment, check if the installed versions of RLiesel and Liesel are compatible with `check_liesel_version()`.
``` r
library(VGAM)
```
Loading required package: stats4
Loading required package: splines
``` r
set.seed(1337)
n <- 1000
x0 <- runif(n)
x1 <- runif(n)
x2 <- runif(n)
y <- rgev(
n,
location = 0 + sin(2 * pi * x0),
scale = exp(-3 + x1),
shape = 0.1 + x2
)
plot(y)
```
``` r
model <- liesel(
response = y,
distribution = "GeneralizedExtremeValue",
predictors = list(
loc = predictor(~ s(x0)),
scale = predictor(~ x1, inverse_link = "Exp"),
concentration = predictor(~ x2)
)
)
```
Installed Liesel version 0.2.4 is compatible.
Now, we can continue in Python and use the `lsl.dist_reg_mcmc()`
function to set up a sampling algorithm with IWLS kernels for the
regression coefficients ($\boldsymbol{\beta}$) and a Gibbs kernel for
the smoothing parameter ($\tau^2$) of the spline. Note that we need to
set $\beta_0$ for $\xi$ to 0.1 manually, because $\xi = 0$ breaks the
sampler.
``` python
import liesel.model as lsl
import jax.numpy as jnp
model = r.model
# concentration == 0.0 seems to break the sampler
model.vars["concentration_p0_beta"].value = jnp.array([0.1, 0.0])
builder = lsl.dist_reg_mcmc(model, seed=42, num_chains=4)
builder.set_duration(warmup_duration=1000, posterior_duration=1000)
engine = builder.build()
```
liesel.goose.engine - INFO - Initializing kernels...
liesel.goose.engine - INFO - Done
``` python
engine.sample_all_epochs()
```
liesel.goose.engine - INFO - Starting epoch: FAST_ADAPTATION, 75 transitions, 25 jitted together
liesel.goose.engine - WARNING - Errors per chain for kernel_00: 7, 12, 10, 4 / 75 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_03: 1, 0, 0, 0 / 75 transitions
liesel.goose.engine - INFO - Finished epoch
liesel.goose.engine - INFO - Starting epoch: SLOW_ADAPTATION, 25 transitions, 25 jitted together
liesel.goose.engine - WARNING - Errors per chain for kernel_00: 2, 1, 1, 3 / 25 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_01: 3, 1, 1, 0 / 25 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_03: 0, 1, 1, 2 / 25 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_04: 2, 1, 0, 1 / 25 transitions
liesel.goose.engine - INFO - Finished epoch
liesel.goose.engine - INFO - Starting epoch: SLOW_ADAPTATION, 50 transitions, 25 jitted together
liesel.goose.engine - WARNING - Errors per chain for kernel_00: 2, 1, 3, 1 / 50 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_01: 1, 1, 1, 2 / 50 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_03: 2, 1, 1, 2 / 50 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_04: 0, 0, 1, 1 / 50 transitions
liesel.goose.engine - INFO - Finished epoch
liesel.goose.engine - INFO - Starting epoch: SLOW_ADAPTATION, 100 transitions, 25 jitted together
liesel.goose.engine - WARNING - Errors per chain for kernel_00: 1, 2, 1, 3 / 100 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_01: 1, 1, 1, 2 / 100 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_03: 1, 0, 1, 2 / 100 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_04: 1, 1, 0, 1 / 100 transitions
liesel.goose.engine - INFO - Finished epoch
liesel.goose.engine - INFO - Starting epoch: SLOW_ADAPTATION, 200 transitions, 25 jitted together
liesel.goose.engine - WARNING - Errors per chain for kernel_00: 1, 1, 1, 1 / 200 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_01: 1, 1, 1, 3 / 200 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_03: 1, 1, 1, 2 / 200 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_04: 3, 0, 1, 1 / 200 transitions
liesel.goose.engine - INFO - Finished epoch
liesel.goose.engine - INFO - Starting epoch: SLOW_ADAPTATION, 500 transitions, 25 jitted together
liesel.goose.engine - WARNING - Errors per chain for kernel_00: 4, 3, 4, 1 / 500 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_01: 1, 2, 1, 1 / 500 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_03: 1, 1, 1, 2 / 500 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_04: 1, 1, 1, 1 / 500 transitions
liesel.goose.engine - INFO - Finished epoch
liesel.goose.engine - INFO - Starting epoch: FAST_ADAPTATION, 50 transitions, 25 jitted together
liesel.goose.engine - WARNING - Errors per chain for kernel_00: 0, 1, 1, 1 / 50 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_01: 1, 2, 1, 0 / 50 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_03: 1, 1, 0, 1 / 50 transitions
liesel.goose.engine - WARNING - Errors per chain for kernel_04: 1, 0, 0, 1 / 50 transitions
liesel.goose.engine - INFO - Finished epoch
liesel.goose.engine - INFO - Finished warmup
liesel.goose.engine - INFO - Starting epoch: POSTERIOR, 1000 transitions, 25 jitted together
liesel.goose.engine - WARNING - Errors per chain for kernel_00: 1, 0, 0, 4 / 1000 transitions
liesel.goose.engine - INFO - Finished epoch
Some tabular summary statistics of the posterior samples:
``` python
import liesel.goose as gs
results = engine.get_results()
gs.Summary(results)
```
Parameter summary:
|
|
|
kernel
|
mean
|
sd
|
q_0.05
|
q_0.5
|
q_0.95
|
sample_size
|
ess_bulk
|
ess_tail
|
rhat
|
|
parameter
|
index
|
|
|
|
|
|
|
|
|
|
|
|
concentration_p0_beta
|
(0,)
|
kernel_00
|
0.068
|
0.051
|
-0.015
|
0.068
|
0.153
|
4000
|
274.887
|
696.928
|
1.012
|
|
(1,)
|
kernel_00
|
1.070
|
0.100
|
0.909
|
1.071
|
1.233
|
4000
|
136.495
|
550.203
|
1.020
|
|
loc_np0_beta
|
(0,)
|
kernel_03
|
0.598
|
0.226
|
0.228
|
0.596
|
0.972
|
4000
|
45.636
|
86.072
|
1.078
|
|
(1,)
|
kernel_03
|
0.322
|
0.124
|
0.115
|
0.324
|
0.521
|
4000
|
96.395
|
249.120
|
1.052
|
|
(2,)
|
kernel_03
|
-0.365
|
0.128
|
-0.614
|
-0.352
|
-0.173
|
4000
|
46.097
|
69.378
|
1.082
|
|
(3,)
|
kernel_03
|
0.362
|
0.065
|
0.252
|
0.363
|
0.471
|
4000
|
62.370
|
152.229
|
1.053
|
|
(4,)
|
kernel_03
|
-0.241
|
0.085
|
-0.381
|
-0.243
|
-0.099
|
4000
|
54.831
|
159.130
|
1.097
|
|
(5,)
|
kernel_03
|
0.170
|
0.031
|
0.119
|
0.168
|
0.221
|
4000
|
35.365
|
145.434
|
1.119
|
|
(6,)
|
kernel_03
|
6.028
|
0.038
|
5.970
|
6.026
|
6.091
|
4000
|
72.399
|
260.400
|
1.032
|
|
(7,)
|
kernel_03
|
0.547
|
0.071
|
0.431
|
0.549
|
0.663
|
4000
|
34.101
|
190.253
|
1.124
|
|
(8,)
|
kernel_03
|
1.703
|
0.030
|
1.656
|
1.702
|
1.751
|
4000
|
69.504
|
213.189
|
1.034
|
|
loc_np0_tau2
|
()
|
kernel_02
|
6.286
|
5.872
|
2.419
|
5.088
|
13.443
|
4000
|
3623.338
|
3754.696
|
1.000
|
|
loc_p0_beta
|
(0,)
|
kernel_04
|
0.027
|
0.003
|
0.023
|
0.027
|
0.032
|
4000
|
76.520
|
49.237
|
1.059
|
|
scale_p0_beta
|
(0,)
|
kernel_01
|
-3.063
|
0.060
|
-3.159
|
-3.066
|
-2.962
|
4000
|
84.772
|
150.819
|
1.045
|
|
(1,)
|
kernel_01
|
1.041
|
0.076
|
0.915
|
1.043
|
1.165
|
4000
|
155.422
|
409.469
|
1.014
|
Error summary:
|
|
|
|
|
count
|
relative
|
|
kernel
|
error_code
|
error_msg
|
phase
|
|
|
|
kernel_00
|
90
|
nan acceptance prob
|
warmup
|
73
|
0.018
|
|
posterior
|
5
|
0.001
|
|
kernel_01
|
90
|
nan acceptance prob
|
warmup
|
30
|
0.008
|
|
posterior
|
0
|
0.000
|
|
kernel_03
|
90
|
nan acceptance prob
|
warmup
|
28
|
0.007
|
|
posterior
|
0
|
0.000
|
|
kernel_04
|
90
|
nan acceptance prob
|
warmup
|
20
|
0.005
|
|
posterior
|
0
|
0.000
|
And the corresponding trace plots:
``` python
fig = gs.plot_trace(results, "loc_p0_beta")
```
/opt/hostedtoolcache/Python/3.10.12/x64/lib/python3.10/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
self._figure.tight_layout(*args, **kwargs)
/opt/hostedtoolcache/Python/3.10.12/x64/lib/python3.10/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
self._figure.tight_layout(*args, **kwargs)
``` python
fig = gs.plot_trace(results, "loc_np0_tau2")
```
/opt/hostedtoolcache/Python/3.10.12/x64/lib/python3.10/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
self._figure.tight_layout(*args, **kwargs)
/opt/hostedtoolcache/Python/3.10.12/x64/lib/python3.10/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
self._figure.tight_layout(*args, **kwargs)
``` python
fig = gs.plot_trace(results, "loc_np0_beta")
```
/opt/hostedtoolcache/Python/3.10.12/x64/lib/python3.10/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
self._figure.tight_layout(*args, **kwargs)
/opt/hostedtoolcache/Python/3.10.12/x64/lib/python3.10/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
self._figure.tight_layout(*args, **kwargs)
``` python
fig = gs.plot_trace(results, "scale_p0_beta")
```
/opt/hostedtoolcache/Python/3.10.12/x64/lib/python3.10/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
self._figure.tight_layout(*args, **kwargs)
/opt/hostedtoolcache/Python/3.10.12/x64/lib/python3.10/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
self._figure.tight_layout(*args, **kwargs)
``` python
fig = gs.plot_trace(results, "concentration_p0_beta")
```
/opt/hostedtoolcache/Python/3.10.12/x64/lib/python3.10/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
self._figure.tight_layout(*args, **kwargs)
/opt/hostedtoolcache/Python/3.10.12/x64/lib/python3.10/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
self._figure.tight_layout(*args, **kwargs)
We need to reset the index of the summary data frame before we can
transfer it to R.
``` python
summary = gs.Summary(results).to_dataframe().reset_index()
```
After transferring the summary data frame to R, we can process it with
packages like dplyr and ggplot2. Here is a visualization of the
estimated spline vs. the true function:
``` r
library(dplyr)
```
Attaching package: 'dplyr'
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
``` r
library(ggplot2)
library(reticulate)
summary <- py$summary
beta <- summary %>%
filter(variable == "loc_np0_beta") %>%
group_by(var_index) %>%
summarize(mean = mean(mean)) %>%
ungroup()
beta <- beta$mean
X <- py_to_r(model$vars["loc_np0_X"]$value)
estimate <- X %*% beta
true <- sin(2 * pi * x0)
ggplot(data.frame(x0 = x0, estimate = estimate, true = true)) +
geom_line(aes(x0, estimate), color = palette()[2]) +
geom_line(aes(x0, true), color = palette()[4]) +
ggtitle("Estimated spline (red) vs. true function (blue)") +
ylab("f") +
theme_minimal()
```
