update examples

Author: Carol-seven

R/gen_prec_sbm.R

@@ -66,7 +66,7 @@
 #' }
 #'
 #' @param weight.paras A list
-#' (default = \code{list(c(shape = 1e4, rate = 1e2), c(min = 0, max = 5))})
+#' (default = \code{list(c(shape = 100, rate = 10), c(min = 0, max = 5))})
 #' specifying the parameters associated with \code{weight.dists}. It must follow
 #' the same length rules as \code{weight.dists}. Each element should be a named
 #' vector or list suitable for the corresponding sampler.
@@ -144,7 +144,7 @@ gen_prec_sbm <- function(d,
                          prob.mat = NULL, within.prob = 0.25, between.prob = 0.05,
                          weight.mat = NULL,
                          weight.dists = list("gamma", "unif"),
-                         weight.paras = list(c(shape = 1e4, rate = 1e2),
+                         weight.paras = list(c(shape = 100, rate = 10),
                                              c(min = 0, max = 5)),
                          cond.target = 100) {
 

R/performance.R

@@ -11,7 +11,8 @@
 #' (typically true) precision matrix.
 #'
 #' @return
-#' A data frame with one row per performance metric and two columns:
+#' A data frame of S3 class \code{"performance"}, with one row per performance
+#' metric and two columns:
 #' \describe{
 #' \item{measure}{The name of each performance metric. The reported metrics
 #' include: sparsity, Frobenius norm loss, Kullback-Leibler divergence,
@@ -89,6 +90,9 @@
 #' \tab Negative predictive value (NPV) = TN / (TN + FN) = 1 - FOR \tab \tab \cr
 #' }
 #'
+#' @example
+#' inst/example/ex-performance.R
+#'
 #' @export
 
 performance <- function(hatOmega, Omega) {
@@ -143,7 +147,7 @@ performance <- function(hatOmega, Omega) {
   ## FPR <- FP / sum(Omega_edge == 0)
   F1 <- 2 * TP / (2*TP + FN + FP)
   ## F1 <- 2 * precision * recall / (precision + recall)
-  MCC <- (TP * TN - FP * FN) / sqrt((TP + FP) * (TP + FN) * (TN + FP) * (TN + FN))
+  MCC <- (TP * TN - FP * FN) / (sqrt(TP + FP) * sqrt(TP + FN) * sqrt(TN + FP) * sqrt(TN + FN))
   ## MCC <- (TP * TN - FP * FN) / sqrt(
   ##   sum(hatOmega_edge != 0) * sum(Omega_edge != 0) * sum(Omega_edge == 0) * sum(hatOmega_edge == 0)
   ## )
@@ -154,5 +158,14 @@ performance <- function(hatOmega, Omega) {
     value = c(sparsity, FL, KL, QL, SL,
               TP, TN, FP, FN, TPR, FPR, F1, MCC)
   )
+  class(result) <- c("performance", class(result))
   return(result)
 }
+
+
+#' @noRd
+
+print.performance <- function(x, ...) {
+  x$value <- round(x$value, 4)
+  NextMethod("print", x)
+}

inst/example/ex-gen_prec_sbm.R

@@ -6,20 +6,20 @@ set.seed(1234)
 ## block-structured precision matrix based on SBM
 #### case 1: base R distribution
 sim1 <- gen_prec_sbm(d = 100, K = 5,
-                     within.prob = 0.5, between.prob = 0.05,
+                     within.prob = 0.25, between.prob = 0.1,
                      weight.dists = list("gamma", "unif"),
-                     weight.paras = list(c(shape = 20, scale = 5),
-                                         c(min = 0, max = 1)),
+                     weight.paras = list(c(shape = 100, scale = 1e2),
+                                         c(min = 0, max = 10)),
                      cond.target = 100)
 #### visualization
 plot(sim1)
 
 #### case 2: user-defined sampler
 my_gamma <- function(n) {
-  rgamma(n, shape = 10, scale = 5)
+  rgamma(n, shape = 1e4, scale = 1e2)
 }
 sim2 <- gen_prec_sbm(d = 100, K = 5,
-                     within.prob = 0.5, between.prob = 0.05,
+                     within.prob = 0.2, between.prob = 0.05,
                      weight.dists = list(my_gamma, "unif"),
                      weight.paras = list(NULL,
                                          c(min = 0, max = 1)),

inst/example/ex-grasps.R

@@ -4,20 +4,24 @@ library(grasps)
 set.seed(1234)
 
 ## block-structured precision matrix based on SBM
-sim <- gen_prec_sbm(d = 60, K = 3,
-                    within.prob = 0.5, between.prob = 0.05,
+sim <- gen_prec_sbm(d = 30, K = 3,
+                    within.prob = 0.25, between.prob = 0.05,
                     weight.dists = list("gamma", "unif"),
-                    weight.paras = list(c(shape = 20, scale = 5), c(min = 0, max = 1)),
+                    weight.paras = list(c(shape = 100, rate = 10),
+                                        c(min = 0, max = 5)),
                     cond.target = 100)
 ## visualization
 plot(sim)
 
 ## n-by-d data matrix
 library(MASS)
-X <- mvrnorm(n = 30, mu = rep(0, 60), Sigma = sim$Sigma)
+X <- mvrnorm(n = 20, mu = rep(0, 30), Sigma = sim$Sigma)
 
-## lasso, BIC
-res <- grasps(X = X, membership = sim$membership, penalty = "lasso", crit = "BIC")
+## adapt, BIC
+res <- grasps(X = X, membership = sim$membership, penalty = "adapt", crit = "BIC")
 
 ## visualization
 plot(res)
+
+## performance
+performance(hatOmega = res$hatOmega, Omega = sim$Omega)

inst/example/ex-performance.R

@@ -0,0 +1,27 @@
+library(grasps)
+
+## reproducibility for everything
+set.seed(1234)
+
+## block-structured precision matrix based on SBM
+sim <- gen_prec_sbm(d = 30, K = 3,
+                    within.prob = 0.25, between.prob = 0.05,
+                    weight.dists = list("gamma", "unif"),
+                    weight.paras = list(c(shape = 10, scale = 5),
+                                        c(min = 0, max = 5)),
+                    cond.target = 100)
+## visualization
+plot(sim)
+
+## n-by-d data matrix
+library(MASS)
+X <- mvrnorm(n = 20, mu = rep(0, 30), Sigma = sim$Sigma)
+
+## adapt, BIC
+res <- grasps(X = X, membership = sim$membership, penalty = "adapt", crit = "BIC")
+
+## visualization
+plot(res)
+
+## performance
+performance(hatOmega = res$hatOmega, Omega = sim$Omega)

inst/example/ex-plot.blkmat.R

@@ -4,10 +4,11 @@ library(grasps)
 set.seed(1234)
 
 ## block-structured precision matrix based on SBM
-sim <- gen_prec_sbm(d = 100, K = 5,
-                    within.prob = 0.5, between.prob = 0.05,
+sim <- gen_prec_sbm(d = 100, K = 10,
+                    within.prob = 0.2, between.prob = 0.05,
                     weight.dists = list("gamma", "unif"),
-                    weight.paras = list(c(shape = 20, scale = 5), c(min = 0, max = 1)),
+                    weight.paras = list(c(shape = 100, scale = 10),
+                                        c(min = 0, max = 5)),
                     cond.target = 100)
 
 ## visualization