get_elts: The function wrapper to get the elements necessary for...

In genscore: Generalized Score Matching Estimators

The function wrapper to get the elements necessary for calculations for all settings.

Description

The function wrapper to get the elements necessary for calculations for all settings.

Usage

get_elts(
  h_hp,
  x,
  setting,
  domain,
  centered = TRUE,
  profiled_if_noncenter = TRUE,
  scale = "",
  diagonal_multiplier = 1,
  use_C = TRUE,
  tol = .Machine$double.eps^0.5,
  unif_dist = NULL
)

Arguments

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.
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". If domain$type == "simplex", only "log_log" and "log_log_sum0" are supported, and on the other hand "log_log_sum0" is supported for domain$type == "simplex" only.
domain	A list returned from make_domain() that represents the domain.
centered	A boolean, whether in the centered setting(assume \boldsymbol{\mu}=\boldsymbol{\eta}=0) or not. Default to TRUE.
profiled_if_noncenter	A boolean, whether in the profiled setting (\lambda_{\boldsymbol{\eta}}=0) if non-centered. Parameter ignored if centered=TRUE. Default to TRUE. Can only be FALSE if setting == "log_log_sum0" && centered == FALSE.
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".
diagonal_multiplier	A number >= 1, the diagonal multiplier.
use_C	Optional. A boolean, use C (TRUE) or R (FALSE) functions for computation. Default to TRUE. Ignored if setting == "gaussian" && domain$type == "R".
tol	Optional. A positive number. If setting != "gaussian" || domain$type != "R", function stops if any entry if smaller than -tol, and all entries between -tol and 0 are set to tol, for numerical stability and to avoid violating the assumption that h(\mathbf{x})>0 almost surely.
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.

Details

Computes the \boldsymbol{\Gamma} matrix and the \boldsymbol{g} vector for generalized score matching.

Here, \boldsymbol{\Gamma} is block-diagonal, and in the non-profiled non-centered setting, the j-th block is composed of \boldsymbol{\Gamma}_{\mathbf{KK},j}, \boldsymbol{\Gamma}_{\mathbf{K}\boldsymbol{\eta},j} and its transpose, and finally \boldsymbol{\Gamma}_{\boldsymbol{\eta\eta},j}. In the centered case, only \boldsymbol{\Gamma}_{\mathbf{KK},j} is computed. In the profiled non-centered case,

\boldsymbol{\Gamma}_{j}\equiv\boldsymbol{\Gamma}_{\mathbf{KK},j}-\boldsymbol{\Gamma}_{\mathbf{K}\boldsymbol{\eta},j}\boldsymbol{\Gamma}_{\boldsymbol{\eta}\boldsymbol{\eta},j}^{-1}\boldsymbol{\Gamma}_{\mathbf{K}\boldsymbol{\eta}}^{\top}.

Similarly, in the non-profiled non-centered setting, \boldsymbol{g} can be partitioned p parts, each with a p-vector \boldsymbol{g}_{\mathbf{K},j} and a scalar g_{\boldsymbol{\eta},j}. In the centered setting, only \boldsymbol{g}_{\mathbf{K},j} is needed. In the profiled non-centered case,

\boldsymbol{g}_j\equiv\boldsymbol{g}_{\mathbf{K},j}-\boldsymbol{\Gamma}_{\mathbf{K}\boldsymbol{\eta},j}\boldsymbol{\Gamma}_{\boldsymbol{\eta\eta},j}^{-1}g_{\boldsymbol{\eta},j}.

The formulae for the pieces above are

\boldsymbol{\Gamma}_{\mathbf{KK},j}\equiv\frac{1}{n}\sum_{i=1}^nh\left(X_j^{(i)}\right){X_j^{(i)}}^{2a-2}{\boldsymbol{X}^{(i)}}^a{{\boldsymbol{X}^{(i)}}^a}^{\top},

\boldsymbol{\Gamma}_{\mathbf{K}\boldsymbol{\eta},j}\equiv-\frac{1}{n}\sum_{i=1}^nh\left(X_j^{(i)}\right){X_j^{(i)}}^{a+b-2}{\boldsymbol{X}^{(i)}}^a,

\boldsymbol{\Gamma}_{\boldsymbol{\eta\eta},j}\equiv\frac{1}{n}\sum_{i=1}^nh\left(X_j^{(i)}\right){X_j^{(i)}}^{2b-2},

\boldsymbol{g}_{\mathbf{K},j}\equiv\frac{1}{n}\sum_{i=1}^n\left(h'\left(X_j^{(i)}\right){X_j^{(i)}}^{a-1}+(a-1)h\left(X_j^{(i)}\right){X_j^{(i)}}^{a-2}\right){\boldsymbol{X}^{(i)}}^a+ah\left(X_j^{(i)}\right){X_j^{(i)}}^{2a-2}\boldsymbol{e}_{j,p},

\boldsymbol{g}_{\boldsymbol{\eta},j}\equiv\frac{1}{n}\sum_{i=1}^n-h'\left(X_j^{(i)}\right){X_j^{(i)}}^{b-1}-(b-1)h\left(X_j^{(i)}\right){X_j^{(i)}}^{b-2},

where \boldsymbol{e}_{j,p} is the p-vector with 1 at the j-th position and 0 elsewhere.

In the profiled non-centered setting, the function also returns t_1 and t_2 defined as

\boldsymbol{t}_1\equiv\boldsymbol{\Gamma}_{\boldsymbol{\eta\eta}}^{-1}\boldsymbol{g}_{\boldsymbol{\eta}},\quad\boldsymbol{t}_2\equiv\boldsymbol{\Gamma}_{\boldsymbol{\eta\eta}}^{-1}\boldsymbol{\Gamma}_{\mathbf{K}\boldsymbol{\eta}}^{\top},

so that \hat{\boldsymbol{\eta}}=\boldsymbol{t}_1-\boldsymbol{t}_2\mathrm{vec}(\hat{\mathbf{K}}).

Value

A list that contains the elements necessary for estimation.

n	The sample size.
p	The dimension.
centered	The centered setting or not. Same as input.
scale	The scaling method. Same as input.
diagonal_multiplier	The diagonal multiplier. Same as input.
diagonals_with_multiplier	A vector that contains the diagonal entries of \boldsymbol{\Gamma} after applying the multiplier.
domain_type	The domain type. Same as domain$type in the input.
setting	The setting. Same as input.
g_K	The \boldsymbol{g} vector. In the non-profiled non-centered setting, this is the \boldsymbol{g} sub-vector corresponding to \mathbf{K}. A p^2-vector. Not returned if setting == "gaussian" && domain$type == "R" since it is just diag(p).
Gamma_K	The \boldsymbol{\Gamma} matrix with no diagonal multiplier. In the non-profiled non-centered setting, this is the \boldsymbol{\Gamma} sub-matrix corresponding to \mathbf{K}. A vector of length p^2 if setting == "gaussian" && domain$type == "R" or p^3 otherwise.
g_eta	Returned in the non-profiled non-centered setting. The \boldsymbol{g} sub-vector corresponding to \boldsymbol{\eta}. A p-vector. Not returned if setting == "gaussian" && domain$type == "R" since it is just numeric(p).
Gamma_K_eta	Returned in the non-profiled non-centered setting. The \boldsymbol{\Gamma} sub-matrix corresponding to interaction between \mathbf{K} and \boldsymbol{\eta}. If setting == "gaussian" && domain$type == "R", returns a vector of length p, or p^2 otherwise.
Gamma_eta	Returned in the non-profiled non-centered setting. The \boldsymbol{\Gamma} sub-matrix corresponding to \boldsymbol{\eta}. A p-vector. Not returned if setting == "gaussian" && domain$type == "R" since it is just rep(1,p).
t1, t2	Returned in the profiled non-centered setting, where the \boldsymbol{\eta} estimate can be retrieved from \boldsymbol{t_1}-\boldsymbol{t_2}\hat{\mathbf{K}} after appropriate resizing.

If domain$type == "simplex", the following are also returned.

Gamma_K_jp	A matrix of size p by p(p-1). The (j-1)*p+1 through j*p columns represent the interaction matrix between the j-th column and the m-th column of K.
Gamma_Kj_etap	Non-centered only. A matrix of size p by p(p-1). The j-th column represents the interaction between the j-th column of K and eta[p].
Gamma_Kp_etaj	Non-centered only. A matrix of size p by p(p-1). The j-th column represents the interaction between the p-th column of K and eta[j]. Note that it is equal to Gamma_Kj_etap if setting does not contain substring "sum0".
Gamma_eta_jp	Non-centered only. A vector of size p-1. The j-th component represents the interaction between eta[j] and eta[p].

Examples

n <- 30
p <- 10
eta <- rep(0, p)
K <- diag(p)
dm <- 1 + (1-1/(1+4*exp(1)*max(6*log(p)/n, sqrt(6*log(p)/n))))

# Gaussian on R^p:
domain <- make_domain("R", p=p)
x <- mvtnorm::rmvnorm(n, mean=solve(K, eta), sigma=solve(K))
# Equivalently:

x2 <- gen(n, setting="gaussian", abs=FALSE, eta=eta, K=K, domain=domain, finite_infinity=100, 
        xinit=NULL, burn_in=1000, thinning=100, verbose=FALSE)

elts <- get_elts(NULL, x, "gaussian", domain, centered=TRUE, scale="norm", diag=dm)
elts <- get_elts(NULL, x, "gaussian", domain, FALSE, profiled=FALSE, scale="sd", diag=dm)

# Gaussian on R_+^p:
domain <- make_domain("R+", p=p)
x <- tmvtnorm::rtmvnorm(n, mean = solve(K, eta), sigma = solve(K),
       lower = rep(0, p), upper = rep(Inf, p), algorithm = "gibbs",
       burn.in.samples = 100, thinning = 10)
# Equivalently:

x2 <- gen(n, setting="gaussian", abs=FALSE, eta=eta, K=K, domain=domain,
       finite_infinity=100, xinit=NULL, burn_in=1000, thinning=100, verbose=FALSE)

h_hp <- get_h_hp("min_pow", 1, 3)
elts <- get_elts(h_hp, x, "gaussian", domain, centered=TRUE, scale="norm", diag=dm)

# Gaussian on sum(x^2) > 1 && sum(x^(1/3)) > 1 with x allowed to be negative
domain <- make_domain("polynomial", p=p, rule="1 && 2",
       ineqs=list(list("expression"="sum(x^2)>1", abs=FALSE, nonnegative=FALSE),
                      list("expression"="sum(x^(1/3))>1", abs=FALSE, nonnegative=FALSE)))
xinit <- rep(sqrt(2/p), p)
x <- gen(n, setting="gaussian", abs=FALSE, eta=eta, K=K, domain=domain, finite_infinity=100, 
       xinit=xinit, seed=2, burn_in=1000, thinning=100, verbose=FALSE)
h_hp <- get_h_hp("min_pow", 1, 3)
elts <- get_elts(h_hp, x, "gaussian", domain, centered=FALSE,
       profiled_if_noncenter=TRUE, scale="", diag=dm)

# exp on ([0, 1] v [2,3])^p
domain <- make_domain("uniform", p=p, lefts=c(0,2), rights=c(1,3))
x <- gen(n, setting="exp", abs=FALSE, eta=eta, K=K, domain=domain,
       xinit=NULL, seed=2, burn_in=1000, thinning=100, verbose=FALSE)
h_hp <- get_h_hp("min_pow", 1.5, 3)
elts <- get_elts(h_hp, x, "exp", domain, centered=TRUE, scale="", diag=dm)
elts <- get_elts(h_hp, x, "exp", domain, centered=FALSE,
       profiled_if_noncenter=FALSE, scale="", diag=dm)

# gamma on {x1 > 1 && log(1.3) < x2 < 1 && x3 > log(1.3) && ... && xp > log(1.3)}
domain <- make_domain("polynomial", p=p, rule="1 && 2 && 3",
       ineqs=list(list("expression"="x1>1", abs=FALSE, nonnegative=TRUE),
                      list("expression"="x2<1", abs=FALSE, nonnegative=TRUE),
                      list("expression"="exp(x)>1.3", abs=FALSE, nonnegative=TRUE)))
set.seed(1)
xinit <- c(1.5, 0.5, abs(stats::rnorm(p-2))+log(1.3))
x <- gen(n, setting="gamma", abs=FALSE, eta=eta, K=K, domain=domain, finite_infinity=100, 
       xinit=xinit, seed=2, burn_in=1000, thinning=100, verbose=FALSE)
h_hp <- get_h_hp("min_pow", 1.5, 3)
elts <- get_elts(h_hp, x, "gamma", domain, centered=TRUE, scale="", diag=dm)
elts <- get_elts(h_hp, x, "gamma", domain, centered=FALSE,
       profiled_if_noncenter=FALSE, scale="", diag=dm)

# a0.6_b0.7 on {x in R_+^p: sum(log(x))<2 || (x1^(2/3)-1.3x2^(-3)<1 && exp(x1)+2.3*x2>2)}
domain <- make_domain("polynomial", p=p, rule="1 || (2 && 3)",
       ineqs=list(list("expression"="sum(log(x))<2", abs=FALSE, nonnegative=TRUE),
                      list("expression"="x1^(2/3)-1.3x2^(-3)<1", abs=FALSE, nonnegative=TRUE),
                      list("expression"="exp(x1)+2.3*x2^2>2", abs=FALSE, nonnegative=TRUE)))
xinit <- rep(1, p)
x <- gen(n, setting="ab_3/5_7/10", abs=FALSE, eta=eta, K=K, domain=domain, finite_infinity=100, 
       xinit=xinit, seed=2, burn_in=1000, thinning=100, verbose=FALSE)
h_hp <- get_h_hp("min_pow", 1.4, 3)
elts <- get_elts(h_hp, x, "ab_3/5_7/10", domain, centered=TRUE, scale="", diag=dm)
elts <- get_elts(h_hp, x, "ab_3/5_7/10", domain, centered=FALSE,
       profiled_if_noncenter=TRUE, scale="", diag=dm)

# log_log model on {x in R_+^p: sum_j j * xj <= 1}
domain <- make_domain("polynomial", p=p,
       ineqs=list(list("expression"=paste(paste(sapply(1:p,
                           function(j){paste(j, "x", j, sep="")}), collapse="+"), "<1"),
                     abs=FALSE, nonnegative=TRUE)))
x <- gen(n, setting="log_log", abs=FALSE, eta=eta, K=K, domain=domain,
       finite_infinity=100, xinit=NULL, seed=2, burn_in=1000, thinning=100,
       verbose=FALSE)
h_hp <- get_h_hp("min_pow", 2, 3)
elts <- get_elts(h_hp, x, "log_log", domain, centered=TRUE, scale="", diag=dm)
elts <- get_elts(h_hp, x, "log_log", domain, centered=FALSE,
       profiled_if_noncenter=FALSE, scale="", diag=dm)
# Example of using the uniform distance function to boundary as in Liu (2019)
g0 <- function(x) {
       row_min <- apply(x, 1, min)
       row_which_min <- apply(x, 1, which.min)
       dist_to_sum_boundary <- apply(x, 1, function(xx){
                   (1 - sum(1:p * xx)) / sqrt(p*(p+1)*(2*p+1)/6)})
       grad_sum_boundary <- -(1:p) / sqrt(p*(p+1)*(2*p+1)/6)
       g0 <- pmin(row_min, dist_to_sum_boundary)
       g0d <- t(sapply(1:nrow(x), function(i){
          if (row_min[i] < dist_to_sum_boundary[i]){
             tmp <- numeric(ncol(x)); tmp[row_which_min[i]] <- 1
          } else {tmp <- grad_sum_boundary}
          tmp
       }))
       list("g0"=g0, "g0d"=g0d)
}
elts <- get_elts(NULL, x, "exp", domain, centered=TRUE, profiled_if_noncenter=FALSE,
       scale="", diag=dm, unif_dist=g0)

# log_log_sum0 model on the simplex with K having row and column sums 0 (Aitchison model)
domain <- make_domain("simplex", p=p)
K <- -cov_cons("band", p=p, spars=3, eig=1)
diag(K) <- diag(K) - rowSums(K) # So that rowSums(K) == colSums(K) == 0
eigen(K)$val[(p-1):p] # Make sure K has one 0 and p-1 positive eigenvalues
x <- gen(n, setting="log_log_sum0", abs=FALSE, eta=eta, K=K, domain=domain,
       xinit=NULL, seed=2, burn_in=1000, thinning=100, verbose=FALSE)
h_hp <- get_h_hp("min_pow", 2, 3)
h_hp_dx <- h_of_dist(h_hp, x, domain) # h and h' applied to distance from x to boundary

# Does not assume K has 0 row and column sums
elts_simplex_0 <- get_elts(h_hp, x, "log_log", domain, centered=TRUE, profiled=FALSE,
       scale="", diag=1.5)

# If want K to have row sums and column sums equal to 0 (Aitchison); estimate off-diagonals only
elts_simplex_1 <- get_elts(h_hp, x, "log_log_sum0", domain, centered=FALSE,
       profiled=FALSE, scale="", diag=1.5)
# All entries corresponding to the diagonals of K should be 0:
max(abs(sapply(1:p, function(j){c(elts_simplex_1$Gamma_K[j, (j-1)*p+1:p],
       elts_simplex_1$Gamma_K[, (j-1)*p+j])})))
max(abs(diag(elts_simplex_1$Gamma_K_eta)))
max(abs(diag(matrix(elts_simplex_1$g_K, nrow=p))))

genscore documentation built on May 29, 2024, 9 a.m.