Case study: automatic marginal selection and the full Sklar workflow

1. One column, many candidate shapes

gpum_fit_catalog() takes a numeric column, detects its support and its modality, fits every eligible parametric family with the sampler, and ranks them by predictive comparison. Take the waiting times of the Old Faithful geyser, a positive and visibly non-Gaussian variable.

y <- faithful$waiting
result <- gpum_fit_catalog(y, n_iter = 8000, n_chains = 8, seed = 1,
                           backend = "cpu")
result
#> <gpum_catalog>
#>   observations   : 272
#>   support        : positive
#>   modality       : multimodal (dip p = 0.00181, KDE modes = 2)
#>   selected family: gmm2 (by WAIC)
#>   diagnostic     : R-hat 1.001 (converged); KS distance 0.034
#>   ranking:
#>       family criterion value  delta converged
#>         gmm2      WAIC  2078   0.00      TRUE
#>      weibull      WAIC  2174  95.99      TRUE
#>       normal      WAIC  2194 116.25      TRUE
#>    student_t      WAIC  2195 117.52      TRUE
#>        gamma      WAIC  2209 131.71      TRUE
#>     logistic      WAIC  2217 139.58      TRUE
#>    lognormal      WAIC  2220 142.39      TRUE
#>      laplace      WAIC  2244 166.03      TRUE
#>  exponential      WAIC  2863 785.42      TRUE

The printed object is the technical diagnostic: the detected support and modality, the selected family, a convergence and goodness-of-fit verdict (R-hat and the Kolmogorov-Smirnov distance of the fitted CDF against the empirical), and the full ranking. The waiting times are bimodal, the dip test rejects unimodality, and the catalogue selects the two-component Gaussian mixture. The mixture is initialised with its components at the lower and upper data quartiles and its means carry a data-anchored prior, so the components start separated on the two modes and the gradient sampler mixes cleanly within each; the draws are then relabeled so the lower-location component is reported first. Every family is fit in an unconstrained parametrisation, so the sampler never proposes an invalid scale or shape.

2. Visual inspection

The plot() method is the visual companion: the fitted density over the histogram (observed against generated), the QQ plot of empirical against fitted quantiles, and the ranking of the candidates.

plot(result)

The mixture density should sit on both humps of the histogram, the QQ points should track the diagonal, and the ranking bar should show the mixture clearly ahead of the unimodal families that cannot represent two modes.

3. The fitted CDF is the transform to the unit square

The winner carries a fitted cumulative distribution, the probability-integral transform. Where the copula vignette used ranks, here the explicit fitted CDF carries the column onto the unit interval, and the two agree closely when the marginal is well chosen.

u_fitted <- marginal_cdf(result, y)
u_rank <- rank(y) / (length(y) + 1)
plot(u_rank, u_fitted, pch = 16, cex = 0.6,
     col = grDevices::adjustcolor("black", 0.5),
     xlab = "rank transform", ylab = "fitted-CDF transform",
     main = "The two probability-integral transforms agree")
abline(0, 1, col = "firebrick", lwd = 1.5, lty = 2)

4. The full Sklar workflow on two columns

With a fitted marginal per column, the dependence is modelled on the fitted-CDF transforms rather than the ranks, so the marginals and the copula together form one explicit model. Take the two variables of the airquality temperature and ozone pair, each given its own automatically selected marginal, then joined by an automatically selected copula.

aq <- airquality[stats::complete.cases(airquality[, c("Temp", "Ozone")]), ]
m_temp <- gpum_fit_catalog(aq$Temp, n_iter = 6000, n_chains = 8, seed = 1,
                           backend = "cpu")
m_ozone <- gpum_fit_catalog(aq$Ozone, n_iter = 6000, n_chains = 8, seed = 1,
                            backend = "cpu")
c(temperature = m_temp$best, ozone = m_ozone$best)
#> temperature       ozone 
#>   "weibull"     "gamma"

# Transform each column by its fitted CDF, then fit the copula on the
# normal scores of those transforms.
clamp <- function(p) pmin(pmax(p, 1e-4), 1 - 1e-4)
u <- clamp(marginal_cdf(m_temp, aq$Temp))
v <- clamp(marginal_cdf(m_ozone, aq$Ozone))
cop <- gpum_copula(data.frame(u = qnorm(u), v = qnorm(v)),
                   family = "auto", n_iter = 8000, n_chains = 8, seed = 1,
                   backend = "cpu")
cop
#> <gpum_copula>
#>   selected family : frank (by WAIC)
#>   theta           : 7.030 [5.481, 8.640]
#>   Kendall's tau   : 0.561 [0.486, 0.625]  (empirical 0.586)
#>   tail dependence : lower 0.000, upper 0.000
#>   diagnostic      : R-hat 1.000 (converged); empirical tau inside the 95% band (adequate fit)
#>   family comparison:
#>    family criterion  value  delta converged
#>     frank      WAIC -95.72  0.000      TRUE
#>  gaussian      WAIC -88.77  6.949      TRUE
#>    gumbel      WAIC -87.56  8.159      TRUE
#>   clayton      WAIC -60.74 34.989      TRUE

The joint model now has three named pieces: the marginal of temperature, the marginal of ozone, and the copula that binds them, each selected by the data. This is the complete factorisation Sklar’s theorem promises, and it is what the synthesis release will sample from to generate synthetic records.

5. Where this stops

The catalogue here fits one column at a time and covers the common continuous families plus the two-component Gaussian mixture; the full catalogue design lists more families and the discrete tier. The mixture handles bimodality; higher-order multimodality needs the three-component mixture, and a dependence structure that is itself a mixture is a copula question, not a marginal one. The fitted CDF is a plug-in at the posterior mean; a fully Bayesian transform would average the CDF over the posterior draws, which the synthesis release will use where the extra fidelity matters.