GEV#

class gpjax.likelihoods.GEV(scale=1.0, shape=0.1, integrator=<gpjax.integrators.GHQuadratureIntegrator object>)[source]#

Bases: AbstractLikelihood

Generalised extreme value (GEV) likelihood for block maxima.

The latent Gaussian process gives the GEV location, \(\mu(x) = f(x)\), and the scale \(\sigma\) and shape \(\xi\) are parameters that all inputs share:

\[y \mid f(x) \sim \operatorname{GEV}(\mu = f(x), \sigma, \xi).\]

This is the latent Gaussian model for spatial extremes (Cooley, Nychka & Naveau, 2007; Davison, Padoan & Ribatet, 2012). The data are block maxima, for example annual maxima of daily rainfall. The shape uses the climate convention: \(\xi > 0\) is a heavy upper tail, \(\xi = 0\) is the Gumbel distribution and \(\xi < 0\) is a bounded upper tail. See GeneralisedExtremeValue.

prior * GEV() gives a NonConjugateModel. Fit it by MAP with log_posterior_density(), or fit a variational family with elbo().

Support. The GEV has a finite end point for \(\xi \neq 0\), and the exact log density is \(-\infty\) outside the support. Thus link_function(), which log_posterior_density, elbo and the Gauss–Hermite expected log-likelihood use, returns the distribution with support_continuation=True: a pseudo-log-density that is finite and has a gradient that points into the support. predict() returns the exact distribution.

Fitting. Prefer a variational family with elbo. When each latent value has one observation, the joint mode of the latent values, the scale and the shape does not exist: the latent can interpolate the data, and then the scale goes to 0 or the shape goes to its upper bound. Thus MAP with log_posterior_density gives useful estimates of the scale and the shape only when inputs repeat, for example many years at each station, where all the years of a station share one latent value.

Starting values. The GEV log density has small gradients far from the data in the heavy tail, so start the latent location and the scale near the data. The Gumbel moment estimates are \(\sigma_0 = s\sqrt{6}/\pi\) and \(\mu_0 = \bar{y} - 0.5772\,\sigma_0\), where \(\bar{y}\) and \(s\) are the sample mean and standard deviation. Use a WhitenedVariationalGaussian, whose latent starts at the mean function, with a Constant mean at \(\mu_0\):

sigma0 = jnp.std(y) * jnp.sqrt(6.0) / jnp.pi
mu0 = jnp.mean(y) - 0.5772 * sigma0
prior = gpx.gps.Prior(
    mean_function=gpx.mean_functions.Constant(mu0), kernel=kernel
)
q = gpx.variational_families.WhitenedVariationalGaussian(
    model=prior * gpx.likelihoods.GEV(scale=sigma0), inducing_inputs=z
)
q = eqx.tree_at(
    lambda q: q.inducing_inputs, q, paramax.non_trainable(q.inducing_inputs)
)

Keep the inducing inputs fixed, as in the last statement. If they train with the step size of the other parameters, the fit can become unstable, and its result can then change with the platform.

The unwhitened VariationalGaussian starts its latent near 0, not at the mean function; give it variational_mean at \(\mu_0\).

Covariates. A covariate in the location, for example \(\mu = \mu_0(x) + \beta\,\mathrm{GMST}(t)\), is an input column: add a Linear kernel with active_dims on that column.

Example

>>> import jax.numpy as jnp
>>> import gpjax as gpx
>>> likelihood = gpx.likelihoods.GEV(scale=12.0, shape=0.1)
>>> prior = gpx.gps.Prior(
...     mean_function=gpx.mean_functions.Constant(40.0),
...     kernel=gpx.kernels.RBF(),
... )
>>> model = prior * likelihood
>>> type(model).__name__
'NonConjugateModel'

Expand for references to gpjax.likelihoods.GEV

GEV

Parameters:

scale (Any)
shape (Any)
integrator (AbstractIntegrator)

link_function(f)[source]#

The GEV distribution of the observations, given the location f.

The distribution uses the support continuation, so its log_prob is finite and has useful gradients for every f; see GeneralisedExtremeValue.

Parameters:: f (Float[Array, ...]) – The latent location values.
Returns:: The distribution of \(y\) given \(f\).
Return type:: GeneralisedExtremeValue

predict(dist)[source]#

The plug-in predictive distribution at the latent mean.

As for Poisson and StudentT, the latent uncertainty is not integrated out: the result is the exact GEV distribution with the location at the latent mean. Use return_level_samples() for return levels that include the latent uncertainty.

Parameters:: dist (MultivariateNormal | GaussianDistribution) – The latent Gaussian process at a finite set of test inputs.
Returns:: The exact GEV, one for each test input.
Return type:: GeneralisedExtremeValue

return_level_samples(dist, return_period, *, key, num_samples=1000)[source]#

Posterior samples of the return level at each test input.

Each sample is \(z_T = f + \sigma c_T(\xi)\), with \(f\) a draw of the latent location from dist and \(c_T\) the standardised return level of gev_return_level(). Take quantiles over the first axis for credible intervals.

The samples include the uncertainty of the latent location. They do not include the uncertainty of the scale and the shape, which are point estimates, or of the kernel hyperparameters. Here the return level is the latent location plus a constant, so its distribution is Gaussian, with the mean and variance of dist shifted by \(\sigma c_T(\xi)\). The samples give the same result, with the same interface as HeteroscedasticGEV.return_level_samples().

If dist has a dense covariance, every sample is one joint draw over the test inputs, so quantities over several inputs (for example a regional mean) are correct. A diagonal covariance (from covariance="diagonal") gives independent draws: the marginal intervals at each input are correct, but a joint quantity is not.

Example

>>> import jax.numpy as jnp
>>> import jax.random as jr
>>> import lineax as lx
>>> import gpjax as gpx
>>> likelihood = gpx.likelihoods.GEV(scale=12.0, shape=0.1)
>>> latent = gpx.distributions.GaussianDistribution(
...     jnp.array([40.0, 45.0]), lx.MatrixLinearOperator(jnp.eye(2))
... )
>>> levels = likelihood.return_level_samples(
...     latent, jnp.array([10.0, 100.0]), key=jr.key(0), num_samples=500
... )
>>> levels.shape
(500, 2, 2)
>>> interval = jnp.quantile(levels, jnp.array([0.05, 0.95]), axis=0)

Parameters:

dist (MultivariateNormal | GaussianDistribution) – The posterior of the latent location at \(N\) test inputs, for example posterior(test_inputs).
return_period (float | int | Float[jaxlib._jax.Array, '...'] | Float[ndarray, '...'] | Int[jaxlib._jax.Array, '...'] | Int[ndarray, '...']) – One return period \(T > 1\), or an array of them, in blocks (years for annual maxima).
key (UInt32[jaxlib._jax.Array, '2'] | Key[jaxlib._jax.Array, '']) – The random key for the latent draws.
num_samples (int) – The number of samples \(S\).

Returns:

The return levels, of shape (S, N) for one return period and (S, *jnp.shape(return_period), N) for an array of them.

Return type:

Float[jaxlib._jax.Array, ‘S …’] | Float[ndarray, ‘S …’]

Expand for references to gpjax.likelihoods.GEV.return_level_samples

GEV

return_period_samples(dist, level, *, key, num_samples=1000)[source]#

Posterior samples of the return period of level at each test input.

Each sample is \(T = 1/P(Y > \text{level})\) for one draw of the latent location; see gev_return_period(). Where the level is above the upper end point of a sample (\(\xi < 0\)), the exceedance probability is 0 and the sample is inf. Keep these samples: jnp.mean(jnp.isinf(samples), axis=0) is the posterior probability that the level is impossible, and 1 / samples gives exceedance probabilities with exact zeros. The uncertainty is as for return_level_samples().

Paired samples. Every sample is one joint draw over the test inputs when dist has a dense covariance. Thus, to compare the same location in two climates, put the location twice in the test inputs, with the two values of the covariate (for example the global mean temperature), and predict with covariance="dense". The ratio of the exceedance probabilities of the two columns is the probability ratio, sample by sample:

x_pair = jnp.array([[site, gmst_past], [site, gmst_now]])
periods = likelihood.return_period_samples(
    posterior(x_pair), level=event_amount, key=key
)
p_past, p_now = 1.0 / periods[:, 0], 1.0 / periods[:, 1]
probability_ratio = p_now / p_past  # inf where p_past == 0

Parameters:

dist (MultivariateNormal | GaussianDistribution) – The posterior of the latent location at \(N\) test inputs.
level (float | int | Float[jaxlib._jax.Array, '...'] | Float[ndarray, '...'] | Int[jaxlib._jax.Array, '...'] | Int[ndarray, '...']) – The level, in the units of the data: a number, or an array that broadcasts to (N,) (for example the observed event amount at each location).
key (UInt32[jaxlib._jax.Array, '2'] | Key[jaxlib._jax.Array, '']) – The random key for the latent draws.
num_samples (int) – The number of samples \(S\).

Returns:

The return periods, of shape (S, N).

Return type:

Float[jaxlib._jax.Array, ‘S N’] | Float[ndarray, ‘S N’]