SLGP with integer outputs

Athénaïs Gautier

Using SLGPs for discrete distributions

In addition to modelling continuous conditional densities, the SLGP implementation can be used to estimate conditional probability mass functions when the response takes values on a discrete support.

For illustration, we return to the quakes catalogue, a dataset shipped with base R (R Core Team 2025): 1000 seismic events of body-wave magnitude above 4.0 recorded since 1964 in the Fiji-Tonga region, originating from the Harvard PRIM-H project. Each event carries its epicentre (lat, long), its hypocentre depth (40-680 km), its magnitude mag, and the number of stations that detected it. Here, we consider earthquake magnitude mag as the response variable. In this dataset, magnitudes are reported on a discretized scale, making them suitable for illustrating the discrete-output formulation.

In this vignette, we demonstrate how to model the conditional distribution

\[ \mathbb{P}(T_x = k \mid X=x), \]

where \(T_x\) denotes the recorded magnitude and \(x\) is the earthquake longitude.

data("quakes")
library(tidyr)
library(dplyr)

# Bin data together for visualisation purpose
df <- quakes %>%
  mutate(long_bin = cut(long, breaks = seq(165, 190, by = 2.5), include.lowest = FALSE)) %>%
  group_by(long_bin) %>%
  mutate(long_bin = paste0(long_bin, "\nn=", n()))%>%
  ungroup()%>%
  mutate(long_bin = factor(long_bin, 
                           levels = sort(unique(long_bin), decreasing = FALSE))) %>%
  data.frame()

range_response <- c(4, 7) 
range_x <- c(165, 190)

We represent the data.

library(ggplot2)
library(ggpubr)
library(viridis)


scatter_plot <- ggplot(df, aes(x = long, y = mag)) +
  geom_point(alpha = 0.5, color = "navy") +
  labs(x = "Longitude (°)",
       y = "Magnitude",
       title = "Observed earthquake magnitudes")  +
  theme_bw()+
  coord_cartesian(xlim=range_x,
                  ylim=range_response)

# Compute normalized frequencies per long_bin
df_bar <- df %>%
  count(long_bin, mag) %>%
  group_by(long_bin) %>%
  mutate(prop = n / sum(n))

# Histogram: Distribution of mag by 'long' bin
hist_plot <- ggplot(df_bar, aes(x = mag)) +
  geom_bar(mapping=aes(y = prop), stat = "identity",
           fill = "darkgrey", color = "grey50", lwd = 0.18, alpha = 0.7)  +
  geom_rug(data = df, aes(x = mag), 
           sides = "b", color = "navy", alpha = 0.5) +
  facet_wrap(~ long_bin, scales = "free_y", nrow=2) +
  labs(x = "Magnitude", 
       y = "Probability density", 
       title = "Histogram of 'magnitude' by 'long' group") +
  theme_bw()+
  coord_cartesian(xlim=range_response,
                  ylim=c(0, 0.5)) 
ggarrange(scatter_plot, hist_plot, ncol = 2, nrow = 1,
          widths = c(0.3, 0.7))
A visual representation of the event magnitudes depending on the longitude in the `quakes` catalogue.

A visual representation of the event magnitudes depending on the longitude in the quakes catalogue.

Fitting a discrete SLGP

The model is fitted in the same way as for a continuous response, except that prediction is performed using the discrete-output option. The response domain is defined by the observed magnitude range, and the SLGP represents a probability mass function over the discrete support.

We use a Random Fourier Feature approximation of a Matérn-\(5/2\) kernel.

Maximum a posteriori estimate

library(SLGP)

modelMAP <- slgp(mag~long, # Use a formula with two indexing variables
                 data=df,
                 method="MAP", #Maximum a posteriori estimation scheme
                 basisFunctionsUsed = "RFF",
                 interpolateBasisFun="WNN", # Accelerate inference
                 hyperparams = list(lengthscale=c(0.1, 0.1), 
                                    sigma2=1), 
                 nIntegral = 31, 
                 sigmaEstimationMethod = "heuristic", 
                 # Set to heuristic for numerical stability                 
                 predictorsLower= c(range_x[1]),
                 predictorsUpper= c(range_x[2]),
                 responseRange= range_response,
                 opts_BasisFun = list(nFreq=150,
                                      MatParam=5/2),
                 discrete=TRUE)

We can represent the conditional densities. We first use the standard plot() method for SLGP objects.

plot( modelMAP,
      newdata = data.frame(long = seq(range_x[1], range_x[2], length.out = 6)),
      draw = "mean",
      panels = TRUE,
      n_response = 31,
      discrete = TRUE)
Conditional magnitude probabilities across longitude under the MAP estimate of the SLGP.

Conditional magnitude probabilities across longitude under the MAP estimate of the SLGP. 


selected_values <- c(167, 180, 185)
gap <- 0.5

df_filtered <- df %>%
  mutate(interval=findInterval(long, c(0, 
                                       selected_values[1]-gap, 
                                       selected_values[1]+gap, 
                                       selected_values[2]-gap, 
                                       selected_values[2]+gap, 
                                       selected_values[3]-gap, 
                                       selected_values[3]+gap)))%>%
  filter(interval %in% c(2, 4, 6))%>%
  group_by(interval)%>%
  mutate(category = paste0("long close to ", c("", selected_values[1],
                                               "", selected_values[2],
                                               "", selected_values[3])[interval], 
                           "\nn=", n()))

names <- sort(unique(df_filtered$category))
dfGrid <- data.frame(expand.grid(selected_values, 
                                 seq(range_response[1], range_response[2],, 31)))
colnames(dfGrid) <- c("long", "mag")
predMAP <- predict(modelMAP, newdata = dfGrid, discrete=TRUE, nIntegral=31)

colnames(predMAP) <- c("long", "mag", "MAP estimator")
predMAP <- predMAP%>%
  pivot_longer(-c("long", "mag"))
predMAP$category <-ifelse(predMAP$long==selected_values[1], names[1],
                          ifelse(predMAP$long==selected_values[2], names[2], names[3]))


df_emp <- df_filtered %>%
  count(category, mag) %>%
  group_by(category) %>%
  mutate(prob = n / sum(n)) %>%
  ungroup()

ggplot(mapping=aes(x = mag)) +
  geom_col(data = df_emp, aes(y = prob), width = 0.09, fill = "darkgrey",
    color = "grey50", linewidth = 0.2, alpha = 0.7)+
  geom_step(data=predMAP, mapping=aes(y=value, group=name, col=name), 
            lwd=1.1, direction = "mid")+
  facet_wrap(~ category, scales = "free_y", nrow=1) +
  labs(x = "Magnitude",
       y = "Probability",    
       title = "Binned 'magnitude' histograms vs SLGP MAP estimates at bins centers") +
  theme_bw()+
  theme(legend.position="bottom",
        legend.direction = "horizontal",
        legend.title = element_blank())+
  coord_cartesian(xlim=range_response,
                  ylim=c(0, 0.2))
Empirical magnitude distributions within longitude bins and SLGP MAP estimates at the corresponding bin centers.

Empirical magnitude distributions within longitude bins and SLGP MAP estimates at the corresponding bin centers.

References

R Core Team. 2025. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. https://www.R-project.org/.