estimate {genscore}

The main function for the generalized score-matching estimator for graphical models.

Description

The main function for the generalized score-matching estimator for graphical models.

Usage

estimate(
  x,
  setting,
  domain,
  elts = NULL,
  centered = TRUE,
  symmetric = "symmetric",
  scale = "",
  lambda1s = NULL,
  lambda_length = NULL,
  lambda_ratio = Inf,
  mode = NULL,
  param1 = NULL,
  param2 = NULL,
  h_hp = NULL,
  unif_dist = NULL,
  verbose = TRUE,
  verbosetext = "",
  tol = 1e-06,
  maxit = 1000,
  BIC_refit = TRUE,
  warmstart = TRUE,
  diagonal_multiplier = NULL,
  eBIC_gammas = c(0, 0.5, 1),
  cv_fold = NULL,
  cv_fold_seed = NULL,
  return_raw = FALSE,
  return_elts = FALSE
)

Arguments

x

An n by p matrix, the data matrix, where n is the sample size and p the dimension.

setting

A string that indicates the distribution type, must be one of "exp", "gamma", "gaussian", "log_log", "log_log_sum0", or of the form "ab_NUM1_NUM2", where NUM1 is the a value and NUM2 is the b value, and NUM1 and NUM2 must be integers or two integers separated by "/", e.g. "ab_2_2", "ab_2_5/4" or "ab_2/3_1/2".

domain

A list returned from make_domain() that represents the domain.

elts

A list (optional), elements necessary for calculations returned by get_elts().

centered

A boolean, whether in the centered setting (assume \boldsymbol{\mu}=\boldsymbol{\eta}=0) or not. Default to TRUE.

symmetric

A string. If equals "symmetric", estimates the minimizer \mathbf{K} over all symmetric matrices; if "and" or "or", use the "and"/"or" rule to get the support. Default to "symmetric".

scale

A string indicating the scaling method. If contains "sd", columns are scaled by standard deviation; if contains "norm", columns are scaled by l2 norm; if contains "center" and setting == "gaussian" && domain$type == "R", columns are centered to have mean zero. Default to "norm".

lambda1s

A vector of lambdas, the penalty parameter for K.

lambda_length

An integer >= 2, the number of lambda1s. Ignored if lambda1s is provided, otherwise a grid of lambdas is automatically chosen so that the results range from an empty graph to a complete graph. Default to 10 if neither lambda1s nor lambda_length is provided.

lambda_ratio

A positive number, the fixed ratio between \lambda_{\mathbf{K}} and \lambda_{\boldsymbol{\eta}}, if \lambda_{\boldsymbol{\eta}}\neq 0 (non-profiled) in the non-centered setting.

mode

A string, the class of the h function. Ignored if elts, or h and hp are provided, or if setting == "gaussian" && domain$type == "R".

param1

A number, the first parameter to the h function. Ignored if elts, or h and hp are provided, or if setting == "gaussian" && domain$type == "R".

param2

A number, the second parameter (may be optional depending on mode) to the h function. Ignored if elts, or h and hp are provided, or if setting == "gaussian" && domain$type == "R".

h_hp

A function that returns a list containing hx=h(x) (element-wise) and hpx=hp(x) (element-wise derivative of h) when applied to a vector or a matrix x, both of which has the same shape as x.

unif_dist

Optional, defaults to NULL. If not NULL, h_hp must be NULL and unif_dist(x) must return a list containing "g0" of length nrow(x) and "g0d" of dimension dim(x), representing the l2 distance and the gradient of the l2 distance to the boundary: the true l2 distance function to the boundary is used for all coordinates in place of h_of_dist; see "Estimating Density Models with Complex Truncation Boundaries" by Liu et al, 2019. That is, (h_j\circ \phi_j)(x_i) in the score-matching loss is replaced by g_0(x_i), the l2 distance of xi to the boundary of the domain.

verbose

Optional. A boolean, whether to output intermediate results.

verbosetext

Optional. A string, text to be added to the end of each printout if verbose == TRUE.

tol

Optional. A number, the tolerance parameter. Default to 1e-6.

maxit

Optional. A positive integer, the maximum number of iterations for each fit. Default to 1000.

BIC_refit

A boolean, whether to get the BIC scores by refitting an unpenalized model restricted to the estimated edges, with lambda1=lambda2=0 and diagonal_multiplier=1. Default to TRUE.

warmstart

Optional. A boolean, whether to use the results from a previous (larger) lambda as a warm start for each new lambda. Default to TRUE.

diagonal_multiplier

A number >= 1, the diagonal multiplier. Optional and ignored if elts is provided. If ncol(x) > ncol(n), a value strictly larger than 1 is recommended. Default to 1+\left(1-\left(1+4e\max\left(6\log p/n, \sqrt{6\log p/n}\right)\right)^{-1}\right).

eBIC_gammas

Optional. A number of a vector of numbers. The \gamma parameter in eBIC. Default to c(0,0.5,1).

cv_fold

Optional. An integer larger than 1 if provided. The number of folds used for cross validation. If provided, losses will be calculated on each fold with model fitted on the other folds, and a lambda_length x cv_fold matrix cv_losses will be returned.

cv_fold_seed

Optional. Seed for generating folds for cross validation.

return_raw

A boolean, whether to return the raw estimates of K. Default to FALSE.

return_elts

A boolean, whether to return the elts used for estimation. Default to FALSE.

Value

edgess

A list of vectors of integers: indices of the non-zero edges.

BICs

A lambda_length by length(eBIC_gammas) matrix of raw eBIC scores (without refitting). If return_raw == FALSE, may contain Infs for rows after the first lambda that gives the complete graph.

lambda1s

A vector of numbers of length lambda_length: the grid of lambda1s over which the estimates are obtained.

converged

A vector of booleans of length lambda_length: indicators of convergence for each fit. If return_raw == FALSE, may contain 0s for all lambdas after the first lambda that gives the complete graph.

iters

A vector of integers of length lambda_length: the number of iterations run for each fit. If return_raw == FALSE, may contain 0s for all lambdas after the first lambda that gives the complete graph.

In addition, if centered == FALSE,

etas

A lambda_length*p matrix of eta estimates with the i-th row corresponding to the i-th lambda1. If return_raw == FALSE, may contain NAs after the first lambda that gives the complete graph.

if centered == FALSE and non-profiled,

lambda2s

A vector of numbers of length lambda_length: the grid of lambda2s over which the estimates are obtained.

if return_raw == TRUE,

raw_estimate

A list that contains lambda_length estimates for K of size ncol(x)*ncol(x).

if BIC_refit == TRUE,

BIC_refits

A lambda_length by length(eBIC_gammas) matrix of refitted eBIC scores, obtained by refitting unpenalized models restricted to the estimated edges. May contain Infs for rows after the first lambda that gives the graph restricted to which an unpenalized model does not have a solution (loss unbounded from below).

if cv_fold is not NULL,

cv_losses

A lambda_length x cv_fold matrix of cross validation losses. If return_raw == FALSE, may contain Infs for all lambdas after the first lambda that gives the complete graph.

if return_elts == TRUE,

elts

A list of elements returned from get_elts().

Examples

# Examples are shown for Gaussian truncated to R+^p only. For other distributions
#   on other types of domains, please refer to \code{gen()} or \code{get_elts()},
#   as the way to call this function (\code{estimate()}) is exactly the same in those cases.
n <- 30
p <- 20
domain <- make_domain("R+", p=p)
mu <- rep(0, p)
K <- diag(p)
lambda1s <- c(0.01,0.1,0.2,0.3,0.4,0.5)
dm <- 1 + (1-1/(1+4*exp(1)*max(6*log(p)/n, sqrt(6*log(p)/n))))
x <- tmvtnorm::rtmvnorm(n, mean = mu, sigma = solve(K),
       lower = rep(0, p), upper = rep(Inf, p), algorithm = "gibbs",
       burn.in.samples = 100, thinning = 10)

## Centered estimates, no elts or h provided, mode and params provided
est1 <- estimate(x, "gaussian", domain=domain, elts=NULL, centered=TRUE,
          symmetric="symmetric", lambda1s=lambda1s, mode="min_pow", 
          param1=1, param2=3, diag=dm, return_raw=TRUE, verbose=FALSE)

h_hp <- get_h_hp("min_pow", 1, 3)
## Centered estimates, no elts provided, h provided; equivalent to est1
est2 <- estimate(x, "gaussian", domain=domain, elts=NULL, centered=TRUE,
          symmetric="symmetric", lambda1s=lambda1s, h_hp=h_hp, diag=dm, 
          return_raw=TRUE, verbose=FALSE)
compare_two_results(est1, est2) ## Should be almost all 0

elts_gauss_c <- get_elts(h_hp, x, setting="gaussian", domain=domain,
            centered=TRUE, diag=dm)
## Centered estimates, elts provided; equivalent to est1 and est2
## Here diagonal_multiplier will be set to the default value, equal to dm above
est3 <- estimate(x, "gaussian", domain=domain, elts=elts_gauss_c,
          symmetric="symmetric", lambda1s=lambda1s, diag=NULL,
          return_raw=TRUE, verbose=FALSE)
compare_two_results(est1, est3) ## Should be almost all 0

## Non-centered estimates with Inf penalty on eta; equivalent to est1~3
est4 <- estimate(x, "gaussian", domain=domain, elts=NULL, centered=FALSE,
          lambda_ratio=0, symmetric="symmetric", lambda1s=lambda1s,
          h=h_hp, diag=dm, return_raw=TRUE, verbose=FALSE)
sum(abs(est4$etas)) ## Should be 0 since non-centered with lambda ratio 0 is equivalent to centered
est4$etas <- NULL ## But different from est1 in that the zero etas are returned in est4
compare_two_results(est1, est4) ## Should be almost all 0


## Profiled estimates, no elts or h provided, mode and params provided
est5 <- estimate(x, "gaussian", domain=domain, elts=NULL, centered=FALSE,
          lambda_ratio=Inf, symmetric="or", lambda1s=lambda1s, mode="min_pow", 
          param1=1, param2=3, diag=dm, return_raw=TRUE, verbose=FALSE)

## Profiled estimates, no elts provided, h provided; equivalent to est5
est6 <- estimate(x, "gaussian", domain=domain, elts=NULL, centered=FALSE,
          lambda_ratio=Inf, symmetric="or", lambda1s=lambda1s,
          h_hp=h_hp, diag=dm, return_raw=TRUE, verbose=FALSE)
compare_two_results(est5, est6) ## Should be almost all 0

elts_gauss_p <- get_elts(h_hp, x, setting="gaussian", domain=domain,
                centered=FALSE, profiled=TRUE, diag=dm)
## Profiled estimates, elts provided; equivalent to est5~6
est7 <- estimate(x, "gaussian", domain=domain, elts=elts_gauss_p, centered=FALSE,
          lambda_ratio=Inf, symmetric="or", lambda1s=lambda1s,
          diagonal_multiplier=NULL, return_raw=TRUE, verbose=FALSE)
compare_two_results(est5, est7) ## Should be almost all 0


## Non-centered estimates, no elts or h provided, mode and params provided
## Using 5-fold cross validation and no BIC refit
est8 <- estimate(x, "gaussian", domain=domain, elts=NULL, centered=FALSE,
          lambda_ratio=2, symmetric="and", lambda_length=100,
          mode="min_pow", param1=1, param2=3, diag=dm, return_raw=TRUE,
          BIC_refit=FALSE, cv_fold=5, cv_fold_seed=2, verbose=FALSE)

## Non-centered estimates, no elts provided, h provided; equivalent to est5
## Using 5-fold cross validation and no BIC refit
est9 <- estimate(x, "gaussian", domain=domain, elts=NULL, centered=FALSE,
          lambda_ratio=2, symmetric="and", lambda_length=100, h_hp=h_hp, diag=dm, 
          return_raw=TRUE, BIC_refit=FALSE, cv_fold=5, cv_fold_seed=2, verbose=FALSE)
compare_two_results(est8, est9) ## Should be almost all 0

elts_gauss_np <- get_elts(h_hp, x, setting="gaussian", domain=domain, centered=FALSE,
                profiled=FALSE, diag=dm)
## Non-centered estimates, elts provided; equivalent to est8~9
## Using 5-fold cross validation and no BIC refit
est10 <- estimate(x, "gaussian", domain, elts=elts_gauss_np, centered=FALSE,
           lambda_ratio=2, symmetric="and", lambda_length=100, diag=NULL,
           return_raw=TRUE, BIC_refit=FALSE, cv_fold=5, cv_fold_seed=2, verbose=FALSE)
compare_two_results(est8, est10) ## Should be almost all 0