--- title: "SUN(3,3) maximum likelihood with analytic gradients" author: "Patrick Brown" date: "`r Sys.Date()`" output: pdf_document vignette: > %\VignetteIndexEntry{SUN(3,3) maximum likelihood with analytic gradients} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4 ) library(admvn) ``` This vignette shows how to obtain the MLE of a SUN$(3,3)$ distribution with analytic gradients from `admvn::dsun_fun()`. Orthant probabilities in dimension $\le 3$ use specialized bivariate/trivariate CDF routines (as in `sn`/`mnormt`), and their gradients are analytic. We compare `stats::optim()` and the C++ trust-region driver `sun_mle()`, and check agreement with `sn::dsun` at the fitted value. ## Parameterization A SUN$(3,3)$ observation is $Y = U \mid V > 0$ with $U,V\in\mathbb{R}^3$. The latent structure is defined on the **correlation scale** of $U$, then scaled by $D_u=\operatorname{diag}(\nu)$ with $\nu_j>0$ the marginal standard deviations. Writing $U_\star$ for the standardized latent (unit variances), \[ \begin{bmatrix} U_\star \\ V \end{bmatrix} \sim \mathcal{N}\!\left( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} I & 0 \\ L & I \end{bmatrix} \begin{bmatrix} C_u & 0 \\ 0 & R \end{bmatrix} \begin{bmatrix} I & L^\top \\ 0 & I \end{bmatrix} \right), \qquad U = \xi + D_u U_\star. \] Equivalently, \[ \operatorname{Var} \begin{pmatrix} U \\ V \end{pmatrix} = \begin{bmatrix} D_u C_u D_u & D_u C_u L^\top \\ L C_u D_u & L C_u L^\top + R \end{bmatrix}. \] Excluding $D_u$, the SUN law is therefore determined by a correlation matrix $C_u$ together with loadings $L$ and residual $R$; entries of $L$ in $[-1,1]$ are the natural regime. The free parameter vector has **21** components and is mapped to the `sn` list `dp = list(xi, Omega, Delta, tau, Gamma)` by `make_sun_params()` as follows. **Location.** $\xi = (\xi_1,\xi_2,\xi_3)^\top$. **Marginal scales.** $\nu_j>0$ are the marginal **standard deviations** of the observed-scale covariance $\Omega = D_u C_u D_u$ with $D_u=\operatorname{diag}(\nu)$, so $\Omega_{jj}=\nu_j^2$. ### Correlation matrices Both $C_u$ (correlation of the standardized latent $U_\star$) and the residual correlation $C$ of $V$ use the same unconstrained map: a unit-row Cholesky factor whose free entries are $\tanh$ of unrestricted parameters. For a $3\times 3$ correlation, three free parameters $z=(z_{12},z_{13},z_{23})$ give \[ \begin{aligned} B &= \begin{bmatrix} 1 & 0 & 0 \\ t_{12} & \sqrt{1-t_{12}^2} & 0 \\ t_{13} & t_{23}\sqrt{1-t_{13}^2} & \sqrt{(1-t_{13}^2)(1-t_{23}^2)} \end{bmatrix}, \\ t_{ij} &= \tanh(z_{ij}), \qquad \operatorname{corr} = B B^\top. \end{aligned} \] - **Observed latent correlation** $C_u$: take $z=(\omega_{12},\omega_{13},\omega_{23})$ and $B=B_u$, so $C_u = B_u B_u^\top$. Then $\Omega = D_u C_u D_u$. - **Residual correlation** $C$: take $z=(\gamma_{12},\gamma_{13},\gamma_{23})$ and $B=B_r$, so $C = B_r B_r^\top$. With $M = L C_u L^\top$ and $D = \operatorname{diag}\!\bigl(1 - \operatorname{diag}(M)\bigr)$ (requiring $\operatorname{diag}(M) < 1$), \[ R = D^{1/2} C D^{1/2}, \qquad \Sigma_{vv} = M + R. \] Then $\Sigma_{vv}$ is a correlation matrix ($\operatorname{diag}(\Sigma_{vv}) = 1$) by construction—no post-hoc rescaling of $V$. (Lowercase $\omega_{ij}$ and $\gamma_{ij}$ are free coordinates; they are distinct from the `sn` matrices $\Omega$ and $\Gamma$ built below.) ### Skew loadings $L$ multiplies $C_u$ (not $\Omega$). All nine entries are free: \[ L = \begin{bmatrix} L_{11} & L_{12} & L_{13} \\ L_{21} & L_{22} & L_{23} \\ L_{31} & L_{32} & L_{33} \end{bmatrix}. \] ($L = 0$ recovers a Gaussian; a small diagonal $L$ is a mild skew start. Default box constraints keep diagonals in $[-1,1]$ and off-diagonals $L_{ij}$, $L_{ji}$ in $[-0.5,0.5]$.) ### `sn` cross-covariance quantities Writing $\Sigma_{uv} = D_u C_u L^\top$, \[ \Delta = D_u^{-1}\Sigma_{uv} = C_u L^\top, \qquad \Gamma = \Sigma_{vv}, \qquad \tau = 0. \] ### Parameter order Names used in code: \[ \begin{aligned} &\xi_1,\xi_2,\xi_3,\; \nu_1,\nu_2,\nu_3,\; \omega_{12},\omega_{13},\omega_{23},\\ &L_{11},L_{12},L_{13},L_{21},L_{22},L_{23},L_{31},L_{32},L_{33},\; \gamma_{12},\gamma_{13},\gamma_{23}. \end{aligned} \] ### Free vector to `dp` (summary) | block | map | |-------|-----| | $\xi$ | identity (location) | | $\nu$ | positive marginal SDs; $\Omega = D_u C_u D_u$ with $D_u=\operatorname{diag}(\nu)$ | | $\omega_{ij}$ | $t=\tanh(\omega)$; unit-row Cholesky $\to C_u$ | | $L_{ij}$ | free loadings on the correlation scale | | $\gamma_{ij}$ | $t=\tanh(\gamma)$; unit-row Cholesky $\to$ residual $C$; $R=D^{1/2} C D^{1/2}$ | Then $\Delta = C_u L^\top$, $\Gamma = L C_u L^\top + R$, $\tau=0$. An alternative joint hyperspherical parameterization of the same SUN law is developed in the `sun_mle_h` vignette. ```{r params} true_params <- c( xi1 = 0, xi2 = 0, xi3 = 0, nu1 = 1, nu2 = 1, nu3 = 1, omega12 = 0.2, omega13 = 0.1, omega23 = 0.2, L11 = 0.4, L12 = 0.1, L13 = 0.1, L21 = 0.1, L22 = 0.4, L23 = 0.1, L31 = 0, L32 = 0.1, L33 = 0.35, gamma12 = 0.3, gamma13 = 0.2, gamma23 = 0.3 ) true_sun_params <- make_sun_params(true_params) str(true_sun_params, max.level = 1) # Sigma_vv = Gamma has unit diagonal by construction diag(true_sun_params$Gamma) ``` ## Simulation A modest sample keeps the vignette fast while still illustrating the API. ```{r sim, eval = requireNamespace("sn", quietly = TRUE)} library(sn) set.seed(0) n <- 500 ysim <- sn::rsun(n, dp = true_sun_params) colnames(ysim) <- paste0("Y", 1:3) head(ysim) ``` ## Build a reusable log-likelihood tape `dsun_fun()` fixes the data and builds per-observation shard tapes. For SUN$(3,3)$ the orthant CDFs use the specialized $d\le 3$ path (QMC settings are unused for those values). ```{r tape, eval = requireNamespace("sn", quietly = TRUE)} # Unused for d<=3 specialized CDF values; kept for API compatibility. qmc_n_points <- 64L qmc_n_shifts <- 2L build_time <- system.time({ tape <- dsun_fun( ysim, true_params, n_points = qmc_n_points, n_shifts = qmc_n_shifts, n_threads = 2L ) })["elapsed"] check <- tape$eval(true_params, log = TRUE, deriv = 1L) sn_ll <- sum(sn::dsun(ysim, dp = true_sun_params, log = TRUE)) c( sn = sn_ll, admvn = check$value, build_sec = unname(build_time), n_threads = tape$n_threads ) ``` ## MLE with `optim()` and an analytic gradient `tape$optim_fns()` returns cached `fn`/`gr` that share one `deriv = 1` evaluation per new parameter value. ```{r mle-setup, eval = requireNamespace("sn", quietly = TRUE)} start_params <- make_sun33_start_from_normal( mu = colMeans(ysim), Sigma = stats::cov(ysim), skew_strength = 0.1 ) bnd <- sun33_bounds() lower <- bnd$lower upper <- bnd$upper # Keep nu floors a bit above the package default for this demo. lower[4:6] <- pmax(lower[4:6], 1e-3) # Rescale free parameters so L-BFGS-B does not take wild Omega steps. parscale <- pmax(abs(as.numeric(start_params)), 0.1) optim_ctrl <- list( fnscale = -1, maxit = 200, parscale = parscale ) obj <- tape$optim_fns() ``` ```{r mle-admvn, eval = requireNamespace("sn", quietly = TRUE)} fit_ad_time <- system.time({ fit_ad <- optim( start_params, fn = obj$fn, gr = obj$gr, method = "L-BFGS-B", lower = lower, upper = upper, control = optim_ctrl ) }) fit_ad$convergence fit_ad$value unname(fit_ad_time["elapsed"]) ``` ## MLE with C++ `sun_mle()` (trustOptim) `sun_mle()` runs a quasi-Newton trust-region loop entirely in C++ (no R fn/gr round-trips). The `nu` parameters are optimized on the log scale. ```{r mle-trust, eval = requireNamespace("sn", quietly = TRUE) && requireNamespace("trustOptim", quietly = TRUE)} fit_trust_time <- system.time({ fit_trust <- sun_mle( tape, start_params, control = list(maxit = 200L, report.freq = 0L, report.level = 0L) ) }) fit_trust$value fit_trust$iterations unname(fit_trust_time["elapsed"]) ``` ## Comparison with `sn::dsun` Score each fit with both `sn::dsun` and `tape$eval` (no separate value-only `sn` optim — that is much slower and not needed to check agreement). Package tests (`tests/testthat/test-dsun.R`, `test-dsun44.R`, `test-sun-kl.R`) also assert that `dsun` / `dsun_fun` / `evaluate_sun33_slice` match `sn::dsun` on mild parameter values. ```{r compare, eval = requireNamespace("sn", quietly = TRUE)} sn_ll_at <- function(par) { sum(sn::dsun(ysim, dp = make_sun_params(par), log = TRUE)) } ad_ll_at <- function(par) { tape$eval(par, log = TRUE, deriv = 0L)$value } compare <- data.frame( method = "ad_optim", sn_loglik = sn_ll_at(fit_ad$par), admvn_loglik = ad_ll_at(fit_ad$par), iterations = unname(fit_ad$counts[["gradient"]]), elapsed_sec = unname(fit_ad_time["elapsed"]), stringsAsFactors = FALSE ) if (exists("fit_trust")) { compare <- rbind( data.frame( method = "ad_trustOptim", sn_loglik = sn_ll_at(fit_trust$par), admvn_loglik = ad_ll_at(fit_trust$par), iterations = as.numeric(fit_trust$iterations), elapsed_sec = unname(fit_trust_time["elapsed"]), stringsAsFactors = FALSE ), compare ) } rownames(compare) <- NULL compare ``` `sn_loglik` and `admvn_loglik` are \(\sum_i \log f(y_i;\hat\theta)\) at each method's reported parameter estimate. With the specialized $d\le 3$ CDF, the two log-densities agree closely. ## Parameter recovery Compare the true free parameters to the start and the MLE (preferring the trustOptim fit when available). Approximate standard errors use a central-difference diagonal of the observed information ($\operatorname{se}_j \approx 1/\sqrt{-H_{jj}}$). (The AD Hessian pattern currently under-reports curvature for some blocks, so finite differences are used for this table.) ```{r recovery, eval = requireNamespace("sn", quietly = TRUE)} par_hat <- as.numeric(if (exists("fit_trust")) fit_trust$par else fit_ad$par) names(par_hat) <- names(true_params) ll0 <- function(p) tape$eval(p, log = TRUE, deriv = 0L)$value f0 <- ll0(par_hat) se <- vapply(seq_along(par_hat), function(j) { h <- pmax(1e-5, 1e-4 * abs(par_hat[j])) p_up <- p_dn <- par_hat p_up[j] <- par_hat[j] + h p_dn[j] <- par_hat[j] - h if (j %in% 4:6) { # keep nu positive if (p_dn[j] <= 1e-8) p_dn[j] <- max(par_hat[j] * 0.5, 1e-6) h_eff2 <- (p_up[j] - par_hat[j]) * (par_hat[j] - p_dn[j]) # for unequal steps use (f+ - f0)/(x+-x0) - (f0 - f-)/(x0-x-) over mean step g_up <- (ll0(p_up) - f0) / (p_up[j] - par_hat[j]) g_dn <- (f0 - ll0(p_dn)) / (par_hat[j] - p_dn[j]) Hjj <- (g_up - g_dn) / (0.5 * ((p_up[j] - par_hat[j]) + (par_hat[j] - p_dn[j]))) } else { Hjj <- (ll0(p_up) - 2 * f0 + ll0(p_dn)) / h^2 } if (is.finite(Hjj) && Hjj < 0) sqrt(-1 / Hjj) else NA_real_ }, numeric(1)) z_score <- (par_hat - as.numeric(true_params)) / se recovery <- data.frame( true = as.numeric(true_params), start = as.numeric(start_params), mle = as.numeric(par_hat), se = se, z = z_score, row.names = names(true_params) ) knitr::kable(recovery, digits = 3) ``` Free-parameter recovery is often weak: several block-Cholesky coordinates can represent nearly the same SUN law, so individual $L$ / $\gamma$ entries need not match the truth even when the fitted density does (see the log-likelihood comparison above). Location, scales, and the leading $\omega$ coordinates are usually better behaved; wide SEs flag the weakly identified slots. ## Gradient check at the fitted MLE Compare analytic/`CppAD` partials to central finite differences along a few one-dimensional slices through $\hat\theta$. ```{r grad-check, eval = requireNamespace("sn", quietly = TRUE), fig.width = 7, fig.height = 5} par_hat <- as.numeric(if (exists("fit_trust")) fit_trust$par else fit_ad$par) par_names <- names(true_params) if (is.null(par_names)) { par_names <- paste0("p", seq_along(par_hat)) } names(par_hat) <- par_names idx <- match(c("xi1", "nu1", "L11", "gamma12"), par_names) ll <- function(p) tape$eval(p, log = TRUE, deriv = 0L)$value ad_grad_j <- function(p, j) tape$eval(p, log = TRUE, deriv = 1L)$gradient[j] fd_grad_j <- function(p, j, h) { p_up <- p p_dn <- p p_up[j] <- p[j] + h p_dn[j] <- p[j] - h (ll(p_up) - ll(p_dn)) / (2 * h) } op <- par(mfrow = c(2, 2), mar = c(3.2, 3.2, 2.2, 0.8), oma = c(0, 0, 1.5, 0)) on.exit(par(op), add = TRUE) for (j in idx) { center <- par_hat[j] span <- pmax(0.15 * abs(center), 0.05) grid <- center + seq(-span, span, length.out = 9L) if (par_names[j] %in% c("nu1", "nu2", "nu3")) { grid <- grid[grid > 1e-3] } ad_vals <- numeric(length(grid)) fd_vals <- numeric(length(grid)) for (k in seq_along(grid)) { p <- par_hat p[j] <- grid[k] h <- pmax(1e-5, 1e-4 * abs(grid[k])) ad_vals[k] <- ad_grad_j(p, j) fd_vals[k] <- fd_grad_j(p, j, h) } ylim <- range(c(ad_vals, fd_vals), finite = TRUE) plot( grid, ad_vals, type = "l", lwd = 2, col = "#1b4f72", xlab = par_names[j], ylab = expression(partialdiff * ell / partialdiff * theta), main = par_names[j], ylim = ylim ) lines(grid, fd_vals, lwd = 2, col = "#b03a2e", lty = 2) abline(v = center, col = "gray50", lty = 3) } legend( "topleft", legend = c("AD", "finite diff"), col = c("#1b4f72", "#b03a2e"), lty = c(1, 2), lwd = 2, bty = "n", cex = 0.8 ) ```