get_elts: 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
)

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].

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}}).$

Examples

Run this code

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:
# \donttest{
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:
# \donttest{
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))))

Run the code above in your browser using DataLab