# Solve the MDYPL state evolution equations with or without intercept, with signal strength or contaminated signal strength ```r solve_se( kappa, ss, alpha, intercept = NULL, start, corrupted = FALSE, gh = NULL, prox_tol = 1e-10, transform = TRUE, init_method = "Nelder-Mead", init_iter = 50, ... ) ``` ## Arguments - `kappa`: asymptotic ratio of columns/rows of the design matrix. `kappa` should be in `(0, 1)`. - `ss`: square root of signal strength or of corrupted signal strength, depending on whether `corrupted = TRUE` or not. See Details. - `alpha`: the shrinkage parameter of the MDYPL estimator. `alpha` should be in `(0, 1)`. - `intercept`: if `NULL` (default) then the MDYPL state evolution equations for the model with no intercept parameter are solved. If a real then the equations for the models with intercept parameter equal to `intercept` are solved. See Details. - `start`: a vector with starting values for `mu`, `b`,`sigma` (and `iota` if `intercept` is numeric). - `corrupted`: if `FALSE` (default) then `ss` is the square root of the signal strength and `intercept`, if numeric, is the oracle intercept value. If `TRUE`, then `ss` is the square root of the corrupted signal strength, and `intercept`, if numeric, is the limit of the estimator computed by `mdyplFit()` with shrinkage parameter `alpha`. See Details. - `gh`: A list with the Gauss-Hermite quadrature nodes and weights, as returned from `statmod::gauss.quad()` with `kind = "hermite"`. Default is `NULL`, in which case `gh` is set to `statmod::gauss.quad(200, kind = "hermite")`. - `prox_tol`: tolerance for the computation of the proximal operator; default is `1e-10`. - `transform`: if `TRUE` (default), the optimization is with respect to `log(mu)`, `log(b)`,`log(sigma)`, (and `iota` if `intercept` is numeric). If `FALSE`, then it is over `mu`, `b`, `sigma` (and `iota` if `intercept` is numeric). The solution is returned in terms of the latter set, regardless of how optimization took place. - `init_method`: The method to be passed to `optim()`. Default is `"Nelder-Mead"`. - `init_iter`: how many iterations of `optim()` should we make to get starting values for `nleqslv::nleqslv()`? Default is `50`, but can also be `0` in which case `start` is directly passed to `nleqslv:nleqslv()`. `init_iter = "only"` results in only `optim()` being used. See Details. - `...`: further arguments to be passed to `nleqslv::nleqslv()`, unless `init_iter = "only"`, in which case `...` is further arguments to be passed to `optim()`. ## Returns If `intercept = NULL`, a vector with the values of `mu`, `b`,`sigma`. Otherwise, a vector with the values of `mu`, `b`,`sigma`, and `iota`, if `corrupted = FALSE`, or the value of the intercept otherwise. The vector has attributes the state evolution functions at the solution (`"funcs"`), the number of iterations used by the last optimization method (`"iter"`), any messages from the last optimization method (`"message"`), and information on the optimization methods used (`"optimization-chain"`). ## Details `init_iter` iterations of `optim()` with `method = init_method` are used towards minimizing `sum(se)^2`, where `se` is a vector of the state evolution functions. The solution is then passed to `nleqslv::nleqslv()` for a more aggressive iteration. The state evolution equations are given in expressions (8) (model without intercept) and expression (15) (model with intercept) in Sterzinger & Kosmidis (2024). If `corrupted = FALSE` (default), then `ss` is the square root of the signal strength, which is the limit $\gamma^2$ of $var(X \beta)$. If `corrupted = TRUE`, then `ss` is the square root of the corrupted signal strength which is the limit $\nu^2$ of $var(X hat(beta)(\alpha))$, where $hat(\beta)(\alpha)$ is the maximum Diaconis-Ylvisaker prior penalized likelihood (MDYPL) estimator as computed by `mdyplFit()` with shrinkage parameter $alpha$. If `intercept = NULL`, then the state evolution equations are solved for the model without intercept. If `intercept` is a real number, then the state evolution equations for the model with intercept are solved (i.e. with predictor c("$\\eta_i = \\theta +\n$", "$x_i^T \\beta$")). In that case, what `intercept` represents depends on the value of `corrupted`. If `corrupted = FALSE`, `intercept` represents the oracle value of $\theta$, otherwise it represents the limit `iota` of the MDYPL estimator of $\theta$ as computed by `mdyplFit()` with shrinkage parameter `alpha`. Note that `start` is always for `mu`, `b`,`sigma`, as is the result, regardless whether `transform = TRUE` or not. Transformations during optimization are done internally. ## Examples ```r ## Reproducing Table 13 of Zhao et al. (2022, DOI: 10.3150/21-BEJ1401) thetas <- c(0, 0.5, 1, 2, 2.5) gamma0 <- 5 pars3 <- matrix(NA, length(thetas), 3) pars4 <- matrix(NA, length(thetas), 4) colnames(pars4) <- c("I_mu", "I_b", "I_sigma", "I_iota") colnames(pars3) <- c("II_mu", "II_b", "II_sigma") for (i in seq_along(thetas)) { start3 <- c(0.5, 1, 1) pars3[i, ] <- solve_se(kappa = 0.2, ss = sqrt(5 + thetas[i]^2), alpha = 1, start = start3, init_iter = 0) start4 <- c(pars3[i, ], thetas[i]) pars4[i, ] <- solve_se(kappa = 0.2, ss = sqrt(5), intercept = thetas[i], alpha = 1, start = start4, init_iter = 0) } cbind(pars3, pars4) ``` ## References Zhao Q, Sur P, Cand\`es E J (2022). The asymptotic distribution of the MLE in high-dimensional logistic models: Arbitrary covariance. **Bernoulli**, 28 , 1835–1861. tools:::Rd_expr_doi("10.3150/21-BEJ1401") . Sterzinger P, Kosmidis I (2024). Diaconis-Ylvisaker prior penalized likelihood for $p/n \to \kappa \in (0,1)$ logistic regression. **arXiv**:2311.07419v2, [https://arxiv.org/abs/2311.07419](https://arxiv.org/abs/2311.07419). ## Author(s) Ioannis Kosmidis `[aut, cre]` ioannis.kosmidis@warwick.ac.uk , Federico Boiocchi `[ctb]` federico.boiocchi@gmail.com , Philipp Sterzinger `[ctb, earlier Julia code by]` P.Sterzinger@lse.ac.uk

solve_se function