Package 'genscore'

Description

Calculates the area under an ROC curve (AUC).

Usage

AUC(tpfp)

Arguments

Value

A number between 0 and 1, the area under the curve (AUC).

Examples

n <- 40
p <- 50
mu <- rep(0, p)
tol <- 1e-8
K <- cov_cons(mode="sub", p=p, seed=1, spars=0.2, eig=0.1, subgraphs=10)
true_edges <- which(abs(K) > tol & diag(p) == 0)
dm <- 1 + (1-1/(1+4*exp(1)*max(6*log(p)/n, sqrt(6*log(p)/n))))
set.seed(1)
domain <- make_domain("R+", p=p)
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)
est <- estimate(x, setting="gaussian", elts=NULL, domain=domain, centered=TRUE,
         symmetric="symmetric", lambda_length=100, mode="min_pow",
         param1=1, param2=3, diagonal_multiplier=dm)
# Apply tp_fp to each estimated edges set for each lambda
TP_FP <- t(sapply(est$edgess, function(edges){tp_fp(edges, true_edges, p)}))
old.par <- par(mfrow=c(1,1), mar=c(5,5,5,5))
auc <- AUC(TP_FP)
plot(c(), c(),  ylim=c(0,1), xlim=c(0,1), cex.lab=1,
  main=paste("ROC curve, AUC",round(auc,4)), xlab="False Positives",
  ylab="True Positives")
points(TP_FP[,2], TP_FP[,1], type="l")
points(c(0,1), c(0,1), type = "l", lty = 2)
par(old.par)

Description

Takes the vertical average of ROC curves using algorithm 3 from Fawcett (2006). The resulting ROC curve preserves the average AUC.

Usage

avgrocs(rocs, num_true_edges, p)

Arguments

Value

The averaged ROC curve, a matrix with 2 columns and (p^2-p-num_true_edges+1) rows.

References

Fawcett T (2006). “An introduction to ROC analysis.” Pattern Recognition Letters, 27(8), 861–874.

Examples

n <- 40
p <- 50
mu <- rep(0, p)
tol <- 1e-8
domain <- make_domain("R+", p=p)
K <- cov_cons(mode="sub", p=p, seed=1, spars=0.2, eig=0.1, subgraphs=10)
true_edges <- which(abs(K) > tol & diag(p) == 0)
dm <- 1 + (1-1/(1+4*exp(1)*max(6*log(p)/n, sqrt(6*log(p)/n))))
ROCs <- list()
old.par <- par(mfrow=c(2,2), mar=c(5,5,5,5))
for (i in 1:3){
  set.seed(i)
  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)
  est <- estimate(x, setting="gaussian", elts=NULL, domain=domain, centered=TRUE,
           symmetric="symmetric", lambda_length=100, mode="min_pow",
           param1=1, param2=3, diag=dm)
  # Apply tp_fp to each estimated edges set for each lambda
  TP_FP <- t(sapply(est$edgess, function(edges){tp_fp(edges, true_edges, p)}))
  ROCs[[i]] <- TP_FP
  plot(c(), c(),  ylim=c(0,1), xlim=c(0,1), cex.lab=1,
    main=paste("ROC, trial ",i,", AUC ",round(AUC(TP_FP),4),sep=""),
    xlab="False Positives", ylab="True Positives")
  points(TP_FP[,2], TP_FP[,1], type="l")
  points(c(0,1), c(0,1), type = "l", lty = 2)
}
average_ROC <- avgrocs(ROCs, length(true_edges), p)
plot(c(), c(),  ylim=c(0,1), xlim=c(0,1), cex.lab=1,
  main=paste("Average ROC, AUC",round(AUC(average_ROC),4)),
  xlab="False Positives", ylab="True Positives")
points(average_ROC[,2], average_ROC[,1], type="l")
points(c(0,1), c(0,1), type = "l", lty = 2)
par(old.par)

Description

Replaces consecutive "&"s and "|"s in a string to a single "&" and "|".

Usage

beautify_rule(rule)

Arguments

Details

Applied to domain$rule if domain$type == "polynomial".

Value

A string with extra "&"s and "|"s removed.

Examples

beautify_rule("(1 & 2 && 3 &&& 4) | 5 || 6 ||| 7")

Description

Finds the index of the bin a number belongs to using binary search.

Usage

binarySearch_bin(arr, l, r, x)

Arguments

Details

Finds the smallest index i such that arr[i] <= x <= arr[i+1].

Value

The index i such that arr[i] <= x <= arr[i+1].

Examples

binarySearch_bin(1:10, 1, 10, seq(1, 10, by=0.5))
binarySearch_bin(1:10, 5, 8, seq(5, 8, by=0.5))

Description

Calculates penalized or unpenalized loss in K and eta given arbitrary data

Usage

calc_crit(elts, res, penalty)

Arguments

Details

This function calculates the loss in some estimated K and eta given an elts generated using get_elts() with a new dataset x. This is helpful for cross-validation.

Value

A number, the loss.

Examples

# In the following examples, all printed numbers should be close to 0.
# In practice, \code{res} need not be estimates fit to \code{elts},
# but in the examples we use \code{res <- get_results(elts)} just to
# demonstrate that the loss this function returns matches that returned
# by the C code during estimation using \code{get_results}.

n <- 6
p <- 3
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))))

domains <- list(make_domain("R", p=p),
                make_domain("R+", p=p),
                make_domain("uniform", p=p, lefts=c(0,2), rights=c(1,3)),
                make_domain("polynomial", p=p,
                  ineqs=list(list("expression"="sum(x^2)<=1", nonnegative=FALSE, abs=FALSE))))

domains <- c(domains,
             list(make_domain("polynomial", p=p,
                    ineqs=list(list("expression"="sum(x^2)<=1", nonnegative=TRUE, abs=FALSE))),
                  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))),
                  make_domain("simplex", p=p)))

for (domain in domains) {
  if (domain$type == "R" ||
       (domain$type == "uniform" && any(domain$lefts < 0)) ||
       (domain$type == "polynomial" && !domain$ineqs[[1]]$nonnegative))
    settings <- c("gaussian")
  else if (domain$type == "simplex")
    settings <- c("log_log", "log_log_sum0")
  else
    settings <- c("gaussian", "exp", "gamma", "log_log", "ab_3/4_2/3")

  if (domain$type == "simplex")
    symms <- c("symmetric")
  else
    symms <- c("symmetric", "and", "or")

  for (setting in settings) {
    x <- gen(n, setting=setting, 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)

    for (symm in symms) {

       # Centered, penalized loss
       elts <- get_elts(h_hp, x, setting, domain, centered=TRUE, scale="", diag=dm)
       res <- get_results(elts, symm, 0.1)
       print(calc_crit(elts, res, penalty=TRUE) - res$crit) # Close to 0

       # Non-centered, unpenalized loss
       elts_nopen <- get_elts(h_hp, x, setting, domain, centered=TRUE, scale="", diag=1)
       res_nopen <- get_results(elts_nopen, symm, 0)
       print(calc_crit(elts_nopen, res_nopen, penalty=FALSE) - res_nopen$crit) # Close to 0

       # Non-centered, non-profiled, penalized loss
       elts_nc_np <- get_elts(h_hp, x, setting, domain, centered=FALSE,
         profiled_if_noncenter=FALSE, scale="", diag=dm)
       res_nc_np <- get_results(elts_nc_np, symm, lambda1=0.1, lambda2=0.05)
       print(calc_crit(elts_nc_np, res_nc_np, penalty=TRUE) - res_nc_np$crit) # Close to 0

       # Non-centered, non-profiled, unpenalized loss
       elts_nc_np_nopen <- get_elts(h_hp, x, setting, domain, centered=FALSE,
         profiled_if_noncenter=FALSE, scale="", diag=1)
       res_nc_np_nopen <- get_results(elts_nc_np_nopen, symm, lambda1=0, lambda2=0)
       print(calc_crit(elts_nc_np_nopen, res_nc_np_nopen, penalty=FALSE) -
         res_nc_np_nopen$crit) # Close to 0

       if (domain$type != "simplex") {
         # Non-centered, profiled, penalized loss
         elts_nc_p <- get_elts(h_hp, x, setting, domain, centered=FALSE,
           profiled_if_noncenter=TRUE, scale="", diag=dm)
         res_nc_p <- get_results(elts_nc_p, symm, lambda1=0.1)
         if (elts_nc_np$setting != setting || elts_nc_np$domain_type != "R")
           res_nc_p$crit <- res_nc_p$crit - sum(elts_nc_np$g_eta ^ 2 / elts_nc_np$Gamma_eta) / 2
         print(calc_crit(elts_nc_np, res_nc_p, penalty=TRUE) - res_nc_p$crit)  # Close to 0
         # Note that the elts argument cannot be profiled, so
         # calc_crit(elts_nc_p, res_nc_p, penalty=TRUE) is not allowed

         # Non-centered, profiled, unpenalized loss
         elts_nc_p_nopen <- get_elts(h_hp, x, setting, domain, centered=FALSE,
           profiled_if_noncenter=TRUE, scale="", diag=1)
         res_nc_p_nopen <- get_results(elts_nc_p_nopen, symm, lambda1=0)
         if (elts_nc_np_nopen$setting != setting || elts_nc_np_nopen$domain_type != "R")
           res_nc_p_nopen$crit <- (res_nc_p_nopen$crit -
              sum(elts_nc_np_nopen$g_eta ^ 2 / elts_nc_np_nopen$Gamma_eta) / 2)
         print(calc_crit(elts_nc_np_nopen, res_nc_p_nopen, penalty=TRUE) -
           res_nc_p_nopen$crit) # Close to 0
          # Again, calc_crit(elts_nc_p_nopen, res_nc_p, penalty=TRUE) is not allowed
       } # if domain$type != "simplex"

    } # for symm in symms
  } # for setting in settings
} # for domain in domains

Description

Checks if two equally sized numeric vectors satisfy the requirements for being left and right endpoints of a domain defined as a union of intervals.

Usage

check_endpoints(lefts, rights)

Arguments

Details

Both lefts and rights must be non-empty and should have the same length. Suppose lefts and rights both have length l, [lefts[1], rights[1]], ..., [lefts[l], rights[l]] must be an increasing and non-overlapping set of valid intervals, meaning lefts[i] <= rights[i] <= lefts[j] for any i < j (singletons and overlapping at the boundary points are allowed). Inf is not allowed in lefts and -Inf is not allowed in rights.

Value

NULL. Program stops if lefts and rights do not define valid left and right endpoints.

Examples

## [-4,-3], [-2,-1], [0,1], [2,3], [4,5]
check_endpoints(lefts=c(-4,-2,0,2,4), rights=c(-3,-1,1,3,5))
## Not run: 
check_endpoints(lefts=c(), rights=c()) # Cannot be empty
check_endpoints(lefts=c(-4,-2,0,2,4), rights=c(-3,-1,1,3)) # Unequal length
check_endpoints(lefts=c(Inf), rights=c(Inf)) # No Inf in lefts, otherwise invalid interval
check_endpoints(lefts=c(-Inf), rights=c(-Inf)) # No -Inf in rights, otherwise invalid interval
check_endpoints(lefts=c(0, 1), rights=c(2, 3)) # [0,2] and [1,3] overlap, not allowed
check_endpoints(lefts=c(2, 0), rights=c(3, 1)) # [2,3], [0,1] not increasing, not allowed

## End(Not run)
## Singletons and overlapping at the boundary points allowed
check_endpoints(lefts=c(0, 1, 2), rights=c(0, 2, 3))

Description

Compares two lists returned from estimate().

Usage

compare_two_results(res, res2)

Arguments

Value

A list of numbers all of which should be close to 0 if res and res2 are expected to be the same.

Description

Compares two lists returned from get_results().

Usage

compare_two_sub_results(res, res2)

Arguments

Value

A list of numbers all of which should be close to 0 if res and res2 are expected to be the same.

Description

Random generator of inverse covariance matrices.

Usage

cov_cons(mode, p, seed = NULL, spars = 1, eig = 0.1, subgraphs = 1)

Arguments

Details

The function generates an inverse covariance matrix according to the mode argument as follows. The diagonal entries of the matrix are set to the same value such that the minimum eigenvalue of the returned matrix is equal to eig.

"random"

Takes the Q matrix from the QR decomposition of a p by p random matrix with independent $Normal(0,1)$ entries, and calculates $Q' diag(ev) Q$. Randomly zeros out its upper triangular entries using independent uniform Bernoulli(spars) variables, and then symmetrizes the matrix using the upper triangular part.

"sub"

Constructs a block diagonal matrix with subgraphs disconnected subgraphs with equal number of nodes. In each subgraph, takes each entry independently from $Uniform(0.5,1)$, and randomly zeros out its upper triangular entries using independent uniform Bernoulli(spars) variables, and finally symmetrizes the matrix using the upper triangular part. The construction from Section 4.2 of Lin et al. (2016).

"er"

Constructs an Erd\Hos-R\'enyi game with probability spars, and sets the edges to independent $Uniform(0.5,1)$ variables, and finally symmetrizes the matrix using the lower triangular entries.

"band"

Constructs a banded matrix so that the (i,j)-th matrix is nonzero if and only if $|i-j|<=spars$, and is equal to $1-|i-j|/(spars+1)$ if $i!=j$.

"chain"

A chain graph, where the (i,j)-th matrix is nonzero if and only if $|i-j|<=1$, and is equal to 0.5 if $|i-j|==1$. A special case of the "band" construction with spars equal to 1.

Value

A p by p inverse covariance matrix. See details.

References

Lin L, Drton M, Shojaie A (2016). “Estimation of high-dimensional graphical models using regularized score matching.” Electron. J. Stat., 10(1), 806–854.

Examples

p <- 100
K1 <- cov_cons("random", p, seed = 1, spars = 0.05, eig = 0.1)
K2 <- cov_cons("sub", p, seed = 2, spars = 0.5, eig = 0.1, subgraphs=10)
K3 <- cov_cons("er", p, seed = 3, spars = 0.05, eig = 0.1)
K4 <- cov_cons("band", p, spars = 2, eig = 0.1)
K5 <- cov_cons("chain", p, eig = 0.1)

Description

The Cram\'er-Rao lower bound (times n) on the variance for estimating the mean parameter mu from a univariate truncated normal sample, assuming the true variance parameter sigmasq is known.

Usage

crbound_mu(mu, sigmasq)

Arguments

Details

The Cram\'er-Rao lower bound in this case is defined as $\sigma^4/var(X-\mu)$.

Value

A number, the Cram\'er-Rao lower bound.

Description

The Cram\'er-Rao lower bound (times n) on the variance for estimating the variance parameter sigmasq from a univariate truncated normal sample, assuming the true mean parameter mu is known.

Usage

crbound_sigma(mu, sigmasq)

Arguments

Details

The Cram\'er-Rao lower bound in this case is defined as $4\sigma^8/var((X-\mu)^2)$.

Value

A number, the Cram\'er-Rao lower bound .

Description

Computes the sum of absolute differences between two lists using diff_vecs().

Usage

diff_lists(l1, l2, name = NULL)

Arguments

Value

Returns the sum of absolute differences between l1 and l2 if name is NULL, or that between l1[[name]] and l2[[name]] otherwise. If name is not NULL and if name is in exactly one of l1 and l2, returns Inf; if name is in neither, returns NA. Exception: Returns a positive integer if the two elements compared hold NA, NULL or Inf values in different places.

Description

Computes the sum of absolute differences in the finite non-NA/NULL elements between two vectors.

Usage

diff_vecs(l1, l2, relative = FALSE)

Arguments

Value

The sum of (relative) absolute differences in l1 and l2, or a positive integer if two vectors differ in length or hold NA, NULL or Inf values in different places.

Description

Returns a list to be passed to C that represents the domain.

Usage

domain_for_C(domain)

Arguments

Details

Construct a list to be read by C code that represents the domain.

Value

A list of the following elements.

Examples

p <- 30
# The 30-dimensional real space R^30
domain <- make_domain("R", p=p)
domain_for_C(domain)

# The non-negative orthant of the 30-dimensional real space, R+^30
domain <- make_domain("R+", p=p)
domain_for_C(domain)

# x such that sum(x^2) > 10 && sum(x^(1/3)) > 10 with x allowed to be negative
domain <- make_domain("polynomial", p=p, rule="1 && 2",
       ineqs=list(list("expression"="sum(x^2)>10", abs=FALSE, nonnegative=FALSE),
                      list("expression"="sum(x^(1/3))>10", abs=FALSE, nonnegative=FALSE)))
domain_for_C(domain)

# ([0, 1] v [2,3]) ^ p
domain <- make_domain("uniform", p=p, lefts=c(0,2), rights=c(1,3))
domain_for_C(domain)

# x such that {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=FALSE)))
domain_for_C(domain)

# x in R_+^p such that {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)))
domain_for_C(domain)

# x in R_+^p such that {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)))
domain_for_C(domain)

# The (p-1)-simplex
domain <- make_domain("simplex", p=p)
domain_for_C(domain)

# The l-1 ball {sum(|x|) < 1}
domain <- make_domain("polynomial", p=p,
       ineqs=list(list("expression"="sum(x)<1", abs=TRUE, nonnegative=FALSE)))
domain_for_C(domain)

Description

Calculates the eBIC score both with and without refitting an unpenalized model restricted to the estimated support.

Usage

eBIC(res, elts, BIC_refit = TRUE, gammas = c(0, 0.5, 1))

Arguments

Value

A vector of length 2*length(gammas). The first length(gammas) numbers are the eBIC scores without refitting for each gamma value, and the rest are those with refitting if BIC_refit == TRUE, or Inf if BIC_refit == FALSE.

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{eBIC()}) is exactly the same in those cases.
n <- 50
p <- 30
domain <- make_domain("R+", p=p)
h_hp <- get_h_hp("min_pow", 1, 3)
mu <- rep(0, p)
K <- diag(p)
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)

elts_gauss_np <- get_elts(h_hp, x, setting="gaussian", domain=domain,
                centered=FALSE, profiled=FALSE, diag=dm)
res_nc_np <- get_results(elts_gauss_np, symmetric="symmetric",
               lambda1=0.35, lambda2=2, previous_res=NULL,
               is_refit=FALSE)
eBIC(res_nc_np, elts_gauss_np, BIC_refit=TRUE, gammas=c(0,0.5,1))

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

Value

In addition, if centered == FALSE,

if centered == FALSE and non-profiled,

if return_raw == TRUE,

if BIC_refit == TRUE,

if cv_fold is not NULL,

if return_elts == TRUE,

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

Description

Finds the max index in a vector that does not exceed a target number.

Usage

find_max_ind(vals, target, start = 1)

Arguments

Value

The max index i such that vals[i] <= target and i >= start.

Examples

for (i in 1:100) {
   vals <- 1:i
   for (start in 1:i)
      for (target in seq(start, i+0.5, by=0.5))
         if (find_max_ind(vals, target, start) != floor(target))
            stop()
}

Description

Evaluate x^(a/b) and |x|^(a/b) with integer a and b with extension to conventional operations (listed under details) that would otherwise result in NaN.

Usage

frac_pow(x, a, b, abs)

Arguments

Details

Replace x by abs(x) below if abs == TRUE. If a == 0 && b == 0, returns log(x). If a != 0 && b == 0, returns exp(a*x). Otherwise, for b != 0, evaluates x^(a/b) with the following extensions. 0^0 evaluates to 1. If x < 0, returns (-1)^a * |x|^(a/b) if b is odd, or NaN otherwise. If x == 0 && a < 0, returns NaN.

Value

A vector of numbers of the same size as x. See details.

Examples

frac_pow(-5:5, 3, 2, TRUE) - abs(-5:5)^(3/2)
frac_pow(-5:5, 5, 3, FALSE) - sign(-5:5)^5*abs(-5:5)^(5/3)
frac_pow(-5:5, 2, 3, FALSE) - ((-5:5)^2)^(1/3)
frac_pow(c(-5:-1,1:5), 0, 0, TRUE) - log(abs(c(-5:-1,1:5)))
frac_pow(-5:5, 0, 1, FALSE) - 1
frac_pow(-5:5, 3, 0, FALSE) - exp(3*-5:5)

Description

Finds the greatest (positive) common divisor of two integers; if one of them is 0, returns the absolute value of the other number.

Usage

gcd(a, b)

Arguments

Value

The greatest (positive) common divisor of two integers; if one of them is 0, returns the absolute value of the other number.

Examples

gcd(1, 2)
gcd(1, -2)
gcd(12, -18)
gcd(-12, 18)
gcd(15, 0)
gcd(0, -15)
gcd(0, 0)

Description

Random data generator from general a-b graphs with general domain types using adaptive rejection metropolis sampling (ARMS). x^(0/0) treated as log(x) and x^(n/0) as exp(x) for n non-zero. Density only guaranteed to be a proper density when 2*a > b >= 0 or when a = b = 0.

Usage

gen(
  n,
  setting,
  abs,
  eta,
  K,
  domain,
  finite_infinity = NULL,
  xinit = NULL,
  seed = NULL,
  burn_in = 1000,
  thinning = 100,
  verbose = TRUE,
  remove_outofbound = TRUE
)

Arguments

Details

NOTE: For polynomial domains with many ineqs and a rule containing "OR" ("|"), not all samples generated are guaranteed to be inside the domain. It is thus recommended to set remove_outofbound to TRUE and rerun the function with new initial points until the desired number of in-bound samples have been generated.

Randomly generates n samples from the p-variate a-b distributions with parameters $\boldsymbol{\eta}$ and $\mathbf{K}$, where p is the length of $\boldsymbol{\eta}$ or the dimension of the square matrix $\mathbf{K}$.

Letting a=a_numer/a_denom and b=b_numer/b_denom, the a-b distribution is proportional to

$\exp\left(-\frac{1}{2a}{\boldsymbol{x}^a}^{\top}\mathbf{K}{\boldsymbol{x}}^a+\boldsymbol{\eta}^{\top}\frac{\boldsymbol{x}^b-\mathbf{1}_p}{b}\right)$

. Note that x^(0/0) is understood as log(x), and x^(n/0) with nonzero n is exp(n*x), and in both cases the a and b in the denominators in the density are treated as 1.

Value

An $n\times p$ matrix of samples, where $p$ is the length of eta.

Examples

n <- 20
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 sum(x^2) > 10 && sum(x^(1/3)) > 10 with x allowed to be negative
domain <- make_domain("polynomial", p=p, rule="1 && 2",
       ineqs=list(list("expression"="sum(x^2)>10", abs=FALSE, nonnegative=FALSE),
                      list("expression"="sum(x^(1/3))>10", abs=FALSE, nonnegative=FALSE)))
xinit <- rep(sqrt(20/p), p)
x <- gen(n, setting="gaussian", abs=FALSE, eta=eta,  K=K, domain=domain, finite_infinity=100, 
       xinit=xinit, seed=2, burn_in=500, thinning=100, verbose=FALSE)

# 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=500, thinning=100, verbose=TRUE)

# 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=FALSE)))
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=500, thinning=100, verbose=FALSE)

# 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=1e4, 
       xinit=xinit, seed=2, burn_in=500, thinning=100, verbose=FALSE)

# log_log model exp(-log(x) %*% K %*% log(x)/2 + eta %*% log(x)) 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=500, thinning=100, verbose=FALSE)

# log_log 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=500, thinning=100, verbose=FALSE)

# Gumbel_Gumbel model exp(-exp(2x) %*% K %*% exp(2x)/2 + eta %*% exp(-3x)) on {sum(|x|) < 1}
domain <- make_domain("polynomial", p=p,
       ineqs=list(list("expression"="sum(x)<1", abs=TRUE, nonnegative=FALSE)))
K <- diag(p)
x <- gen(n, setting="ab_2/0_-3/0", abs=FALSE, eta=eta, K=K, domain=domain, finite_infinity=100, 
       xinit=NULL, seed=2, burn_in=500, thinning=100, verbose=FALSE)

Description

Analytic solution of the minimized loss for an unpenalized asymmetric model restricted to a given support. Does not work if symmetric == "symmetric".

Usage

get_crit_nopenalty(
  elts,
  exclude = NULL,
  exclude_eta = NULL,
  previous_res = NULL
)

Arguments

Details

If previous_res is provided, exclude and exclude_eta must be NULL or be consistent with the estimated support in previous_res. If previous_res and exclude are both NULL, assume all edges are present. The same applies to the non-profiled non-centered case when previous_res and exclude_eta are both NULL.

Value

A number, the refitted loss.

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{get_crit_nopenalty()}) is exactly the same in those cases.
n <- 50
p <- 30
domain <- make_domain("R+", p=p)
h_hp <- get_h_hp("min_pow", 1, 3)
mu <- rep(0, p)
K <- diag(p)
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)

elts_gauss_np <- get_elts(h_hp, x, setting="gaussian", domain=domain,
                centered=FALSE, profiled=FALSE, diag=dm)
res_nc_np <- get_results(elts_gauss_np, symmetric="symmetric", lambda1=0.35,
               lambda2=2, previous_res=NULL, is_refit=FALSE)
get_crit_nopenalty(elts_gauss_np, previous_res=res_nc_np)

Description

Finds the distance of each element in a matrix x to its boundary of the domain while fixing the others in the same row.

Usage

get_dist(x, domain)

Arguments

Details

Returned matrix dx has its i,j-th component the distance of $x_{i,j}$ to the boundary of domain, assuming $x_{i,-j}$ are fixed. The matrix has the same size of x (n by p), or if domain$type == "simplex" and x has full dimension p, it has p-1 columns. Returned matrix dpx contains the component-wise derivatives of dx in its components. That is, its i,j-th component is 0 if $x_{i,j}$ is unbounded or is bounded from both below and above or is at the boundary, or -1 if $x_{i,j}$ is closer to its lower boundary (or if its bounded from below but unbounded from above), or 1 otherwise.

Value

A list that contains h(dist(x, domain)) and h\'(dist(x, domain)).

Examples

n <- 20
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)

dist <- get_dist(x, domain)
# dx is all Inf and dpx is all 0 since each coordinate is unbounded with domain R
c(all(is.infinite(dist$dx)), all(dist$dpx==0))

# exp on R_+^p:
domain <- make_domain("R+", p=p)
x <- gen(n, setting="exp", abs=FALSE, eta=eta, K=K, domain=domain, finite_infinity=100, 
       xinit=NULL, seed=2, burn_in=1000, thinning=100, verbose=FALSE)
dist <- get_dist(x, domain)
# dx is x and dpx is 1; with domain R+, the distance of x to the boundary is just x itself
c(max(abs(dist$dx - x))<.Machine$double.eps^0.5, all(dist$dpx == 1))

# Gaussian on sum(x^2) > p with x allowed to be negative
domain <- make_domain("polynomial", p=p,
       ineqs=list(list("expression"=paste("sum(x^2)>", p), abs=FALSE, nonnegative=FALSE)))
x <- gen(n, setting="gaussian", abs=FALSE, eta=eta, K=K, domain=domain, finite_infinity=100, 
       xinit=NULL, seed=2, burn_in=1000, thinning=100, verbose=FALSE)
dist <- get_dist(x, domain)
quota <- p - (rowSums(x^2) - x^2) # How much should xij^2 at least be so that sum(xi^2) > p?
# How far is xij from +/-sqrt(quota), if quota >= 0?
dist_to_bound <- abs(x[quota >= 0]) - abs(sqrt(quota[quota >= 0]))
max(abs(dist$dx[is.finite(dist$dx)] - dist_to_bound)) # Should be equal to our own calculations
# dist'(x) should be the same as the sign of x
all(dist$dpx[is.finite(dist$dx)] == sign(x[quota >= 0]))
# quota is negative <-> sum of x_{i,-j}^2 already > p <-> xij unbounded
#   given others <-> distance to boundary is Inf
all(quota[is.infinite(dist$dx)] < 0)

# gamma 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="gamma", abs=FALSE, eta=eta, K=K, domain=domain,
       xinit=NULL, seed=2, burn_in=1000, thinning=100, verbose=FALSE)
dist <- get_dist(x, domain)
# If 0 <= xij <= 1, distance to boundary is min(x-0, 1-x)
max(abs(dist$dx - pmin(x, 1-x))[x >= 0 & x <= 1])
# If 0 <= xij <= 1, dist'(xij) is 1 if it is closer to 0, or -1 if it is closer 1,
#   assuming xij %in% c(0, 0.5, 1) with probability 0
all((dist$dpx == 2 * (1-x > x) - 1)[x >= 0 & x <= 1])
# If 2 <= xij <= 3, distance to boundary is min(x-2, 3-x)
max(abs(dist$dx - pmin(x-2, 3-x))[x >= 2 & x <= 3])
# If 2 <= xij <= 3, dist'(xij) is 1 if it is closer to 2, or -1 if it is closer 3,
#   assuming xij %in% c(2, 2.5, 3) with probability 0
all((dist$dpx == 2 * (3-x > x-2) - 1)[x >= 2 & x <= 3])

# a0.6_b0.7 on {x1 > 1 && 0 < x2 < 1 && x3 > 0 && ... && xp > 0}
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=FALSE)))
set.seed(1)
xinit <- c(1.5, 0.5, abs(stats::rnorm(p-2)) + log(1.3))
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)
dist <- get_dist(x, domain)
# x_{i1} has uniform bound [1, +Inf), so its distance to its boundary is x_{i1} - 1
max(abs(dist$dx[,1] - (x[,1] - 1)))
# x_{i2} has uniform bound [log(1.3), 1], so its distance to its boundary
#   is min(x_{i2} - log(1.3), 1 - x_{i2})
max(abs(dist$dx[,2] - pmin(x[,2] - log(1.3), 1 - x[,2])))
# x_{ij} for j >= 3 has uniform bound [log(1.3), +Inf), so its distance to its boundary
#   is simply x_{ij} - log(1.3)
max(abs(dist$dx[,3:p] - (x[,3:p] - log(1.3))))
# dist\'(xi2) is 1 if it is closer to log(1.3), or -1 if it is closer 1,
#   assuming x_{i2} %in% c(log(1.3), (1+log(1.3))/2, 1) with probability 0
all((dist$dpx[,2] == 2 * (1 - x[,2] > x[,2] - log(1.3)) - 1))
# x_{ij} for j != 2 is bounded from below but unbounded from above, so dist\'(xij) is always 1
all(dist$dpx[,-2] == 1)

# 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)
dist <- get_dist(x, domain)
# Upper bound for j * xij so that sum_j j * xij <= 1
quota <- 1 - (rowSums(t(t(x) * 1:p)) - t(t(x) * 1:p))
# Distance of xij to its boundary is min(xij - 0, quota_{i,j} / j - xij)
max(abs(dist$dx - pmin((t(t(quota) / 1:p) - x), x)))

# log_log 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)
# Note that dist$dx and dist$dpx only has p-1 columns -- excluding the last coordinate in x
dist <- get_dist(x, domain)
# Upper bound for x_{i,j} so that x_{i,1} + ... + x_{i,p-1} <= 1
quota <- 1 - (rowSums(x[,-p]) - x[,-p])
# Distance of x_{i,j} to its boundary is min(xij - 0, quota_{i,j} - xij)
max(abs(dist$dx - pmin(quota - x[,-p], x[,-p])))

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

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.