calc_crit: Calculates penalized or unpenalized loss in K and eta given...

In genscore: Generalized Score Matching Estimators

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

Description

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

Usage

calc_crit(elts, res, penalty)

Arguments

elts	An element list returned from get_elts(). Need not be the same as the elements used to estimate res, but they must be both centered or both non-centered, and their dimension p must match. elts cannot be profiled as this is supposed to be elements for a new data unseen by res, in which case the loss must be explicitly written in K and eta with Gamma and g from a new dataset x.
res	A result list returned from get_results(). Must be centered if elts is centered, and must be non-centered otherwise. Can be profiled. res$p must be equal to elts$p.
penalty	A boolean, indicates whether the loss should be penalized (using elts$diagonals_with_multiplier, res$lambda1 and res$lambda2).

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

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