version 5.0

Author: Marco Scutari

Changelog

@@ -1,6 +1,44 @@
-bnlearn (4.9.4)
-
-  * added some PROTECT()s to pass the CRAN tests.
+bnlearn (5.0)
+
+  * the "effective" argument of nparams() is now deprecated and will be removed
+     by the end of 2025.
+  * the cross-validation loss functions "pred-lw", "cor-lw" and "mse-lw" are
+     now deprecated in favour of "pred", "cor" and "mse" with optional
+     arguments predict = "bayes-lw"; they will be removed in 2025.
+  * the cross-validation loss functions "logl-g" and "logl-cg" are now
+     deprecated in favour of "logl"; they will be removed in 2025.
+  * the cross-validation loss functions "pred", "cor" and "mse" can now be
+     computed with exact inference using predict = "exact".
+  * completed the implementation of KL(), which now supports conditional 
+     Gaussian networks in addition to discrete and Gaussian ones.
+  * implemented Shannon's entropy.
+  * conditional independence tests now have optional arguments like network
+     scores.
+  * added a "custom-test" conditional independence test allowing user-provided
+     test statistics in the same way as "custom" allows user-provided network
+     scores.
+  * the custom score now has label "custom-score", instead of just "custom",
+     for clarity and to make it consistent with "custom-test".
+  * added a "params.threshold" to hard EM methods in bn.fit(), and renamed
+     the log-likelihood threshold to "loglik.threshold".
+  * the log-likelihood stopping rule in hard EM now uses the log-likelihood of
+     the completed data, which works better in the presence of latent
+     variables and is more appropriate accoridng to Koller & Friedman (thanks
+     Laura Azzimonti).
+  * coefficients(), sigma(), fitted() and residuals() return an error when
+     called on bn objects instead of failing silently and returning NULL.
+  * preserve and return the probabilities from predict(..., prob = TRUE) when
+     using the parallel package (thanks Alex Rudge).
+  * logLik() now returns an object of class "logLik" with the expected
+     attributes.
+  * added an identifiable() function to tell whether a bn.fit object contains NA
+     parameter values; and a singular() function to tell whether it is
+     singular (with 0-1 probability distributions, zero standard errors).
+  * the entropy loss in bn.cv() is now estimated now uses the node-average 
+     (log-)likelihood; it does not produce warnings for incomplete data and
+     may occasionally return +Inf instead of NA in some corner cases (for
+     instance, when the model learning from training has NA parameters or
+     conditional probabilities equal to zero).
 
 bnlearn (4.9.3)
 

DESCRIPTION

@@ -2,12 +2,14 @@ Package: bnlearn
 Type: Package
 Title: Bayesian Network Structure Learning, Parameter Learning and
         Inference
-Version: 4.9.4
-Date: 2024-05-02
-Depends: R (>= 4.3.0), methods
+Version: 5.0
+Date: 2024-07-29
+Depends: R (>= 4.4.0), methods
 Suggests: parallel, graph, Rgraphviz, igraph, lattice, gRbase, gRain
-        (>= 1.3-3), ROCR, Rmpfr, gmp
-Author: Marco Scutari [aut, cre], Tomi Silander [ctb], Robert Ness [ctb]
+        (>= 1.3-3), Rmpfr, gmp
+Authors@R: c(person(given = "Marco", family = "Scutari", role = c("aut", "cre"),
+                    email = "[email protected]"),
+             person(given = "Tomi", family = "Silander", role = "ctb"))
 Maintainer: Marco Scutari <[email protected]>
 Description: Bayesian network structure learning, parameter learning and inference.
   This package implements constraint-based (PC, GS, IAMB, Inter-IAMB, Fast-IAMB, MMPC,
@@ -26,6 +28,8 @@ SystemRequirements: USE_C17
 License: GPL (>= 2)
 LazyData: yes
 NeedsCompilation: yes
-Packaged: 2024-05-02 17:13:25 UTC; fizban
+Packaged: 2024-07-30 17:02:42 UTC; fizban
+Author: Marco Scutari [aut, cre],
+  Tomi Silander [ctb]
 Repository: CRAN
-Date/Publication: 2024-05-02 22:30:02 UTC
+Date/Publication: 2024-07-30 17:30:02 UTC

MD5

@@ -26,7 +26,7 @@ export(
   "set.edge", "drop.edge", "add.node", "remove.node", "rename.nodes",
   # get, set and count sets of nodes: parents, children, etc.
   "parents", "parents<-", "children", "children<-", "spouses", "ancestors",
-  "descendants", "root.nodes", "leaf.nodes", "nnodes",
+  "descendants", "root.nodes", "leaf.nodes", "isolated.nodes", "nnodes",
   # get and set adjacency matrices.
   "amat", "amat<-",
   # model string formulas.
@@ -54,19 +54,21 @@ export(
   # data preprocessing and imputation.
   "discretize", "dedup", "impute",
   # inference.
-  "KL", "cpquery",
+  "H", "KL", "cpquery",
   # import/export functions for varous file formats.
   "read.bif", "write.bif", "read.dsc", "write.dsc", "read.net", "write.net",
   "write.dot",
   # utility functions to manipulate test/score counters.
   "test.counter", "increment.test.counter", "reset.test.counter",
   # assorted functions involving network structures.
   "acyclic", "directed", "path.exists", "node.ordering", "subgraph",
+  # assorted functions involving fitted networks.
+  "identifiable", "singular", 
   # assorted functions to extract information.
   "configs", "nparams", "ntests",
   # assorted conversion functions.
   "as.bn", "as.bn.fit", "as.grain", "as.graphNEL", "as.graphAM", "as.igraph",
-  "as.prediction", "as.lm",
+  "as.lm",
   # graph enumeration.
   "count.graphs"
 )
@@ -110,7 +112,6 @@ S3method(as.graphAM, "bn.fit")
 S3method(as.igraph, "bn")
 S3method(as.igraph, "bn.fit")
 S3method(as.bn, "igraph")
-S3method(as.prediction, "bn.strength")
 S3method(as.lm, "bn")
 S3method(as.lm, "bn.fit")
 S3method(as.lm, "bn.fit.gnode")
@@ -120,21 +121,25 @@ S3method(print, "bn.fit.dnode")
 S3method(print, "bn.fit.onode")
 S3method(print, "bn.fit.gnode")
 S3method(print, "bn.fit.cgnode")
+S3method(residuals, "bn")
 S3method(residuals, "bn.fit")
 S3method(residuals, "bn.fit.dnode")
 S3method(residuals, "bn.fit.onode")
 S3method(residuals, "bn.fit.gnode")
 S3method(residuals, "bn.fit.cgnode")
+S3method(fitted, "bn")
 S3method(fitted, "bn.fit")
 S3method(fitted, "bn.fit.dnode")
 S3method(fitted, "bn.fit.onode")
 S3method(fitted, "bn.fit.gnode")
 S3method(fitted, "bn.fit.cgnode")
+S3method(sigma, "bn")
 S3method(sigma, "bn.fit")
 S3method(sigma, "bn.fit.dnode")
 S3method(sigma, "bn.fit.onode")
 S3method(sigma, "bn.fit.gnode")
 S3method(sigma, "bn.fit.cgnode")
+S3method(coef, "bn")
 S3method(coef, "bn.fit")
 S3method(coef, "bn.fit.dnode")
 S3method(coef, "bn.fit.onode")

NAMESPACE

@@ -1,6 +1,7 @@
 
 # compute arcs' strength as the p-value of the test for their removal.
-arc.strength.test = function(network, data, test, alpha, B, debug = FALSE) {
+arc.strength.test = function(network, data, test, alpha, extra.args,
+    debug = FALSE) {
 
   drop = function(arc) {
 
@@ -15,7 +16,7 @@ arc.strength.test = function(network, data, test, alpha, B, debug = FALSE) {
       network$nodes[[arc[2]]]$parents[network$nodes[[arc[2]]]$parents != arc[1]]
 
     a = indep.test(arc[1], arc[2], parents, data = data, test = test,
-          B = B, alpha = alpha)
+          extra.args = extra.args, alpha = alpha)
 
     if (debug) {
 
@@ -228,7 +229,7 @@ arc.strength.boot = function(data, cluster = NULL, R, m, algorithm,
 # compute an approximation of arc and direction strength from the Bayes factors
 # that can be computed from a single MAP network.
 bf.strength.backend = function(x, data, score, extra.args, precBits = 200,
-  debug = FALSE) {
+    debug = FALSE) {
 
   # construct all pairs of nodes.
   nodes = names(x$nodes)

R/arc.strength.R

@@ -2,7 +2,7 @@
 # common steps in constraint-based structure learning: infer arc orientation
 # from the skeleton and the data.
 learn.arc.directions = function(x, cluster = NULL, local.structure, whitelist,
-    blacklist, test, alpha, B = NULL, data, max.sx = ncol(data),
+    blacklist, test, alpha, extra.args = list(), data, max.sx = ncol(data),
     debug = FALSE) {
 
   nodes = names(x)
@@ -18,8 +18,8 @@ learn.arc.directions = function(x, cluster = NULL, local.structure, whitelist,
   # 3.1 detect v-structures.
   vs = do.call("rbind",
          vstruct.detect(nodes = nodes, arcs = arcs, mb = local.structure,
-           data = x, alpha = alpha, B = B, test = test, blacklist = blacklist,
-           max.sx = max.sx, debug = debug))
+           data = x, alpha = alpha, extra.args = extra.args, test = test,
+           blacklist = blacklist, max.sx = max.sx, debug = debug))
   rownames(vs) = NULL
 
   if (!is.null(vs)) {
@@ -31,21 +31,14 @@ learn.arc.directions = function(x, cluster = NULL, local.structure, whitelist,
 
   }#THEN
 
-  # save the status of the learning algorithm.
-  learning = list(whitelist = whitelist, blacklist = blacklist,
-    test = test, args = list(alpha = alpha), ntests = test.counter())
-  # include also the number of permutations/bootstrap samples if it makes sense.
-  if (!is.null(B))
-    learning$args$B = B
-
   # 4. propagate directions.
-  pdag = list(learning = learning,
+  pdag = list(learning = list(whitelist = whitelist, blacklist = blacklist),
               nodes = structure(rep(0, length(nodes)), names = nodes),
               arcs = arcs)
 
-  return(cpdag.backend(pdag, fix = TRUE, debug = debug))
+  return(cpdag.backend(pdag, fix = TRUE, debug = debug)$arcs)
 
-}#SECOND.PRINCIPLE
+}#LEARN.ARC.DIRECTIONS
 
 # construct a fake markov blanket using all the nodes within distance 2.
 fake.markov.blanket = function(learn, target) {
@@ -59,8 +52,9 @@ fake.markov.blanket = function(learn, target) {
 }#FAKE.MARKOV.BLANKET
 
 # build the neighbourhood of a node from the markov blanket.
-neighbour = function(x, mb, data, alpha, B = NULL, whitelist, blacklist,
-  test, empty.dsep = TRUE, markov = TRUE, max.sx = ncol(x), debug = FALSE) {
+neighbour = function(x, mb, data, alpha, extra.args = list(),
+    whitelist, blacklist, test, empty.dsep = TRUE, markov = TRUE,
+    max.sx = ncol(x), debug = FALSE) {
 
   # initialize the neighbourhood using the markov blanket.
   candidate.neighbours = mb[[x]]
@@ -121,7 +115,8 @@ neighbour = function(x, mb, data, alpha, B = NULL, whitelist, blacklist,
 
     a = allsubs.test(x = x, y = y, sx = dsep.set,
           min = ifelse(empty.dsep, 0, 1), max = min(length(dsep.set), max.sx),
-          data = data, test = test, alpha = alpha, B = B, debug = debug)
+          data = data, test = test, alpha = alpha, extra.args = extra.args,
+          debug = debug)
 
     # update the neighbourhood.
     if (a["p.value"] > alpha)
@@ -138,8 +133,8 @@ neighbour = function(x, mb, data, alpha, B = NULL, whitelist, blacklist,
 }#NEIGHBOUR
 
 # detect v-structures in the graph.
-vstruct.detect = function(nodes, arcs, mb, data, alpha, B = NULL, test,
-    blacklist, max.sx = ncol(data), debug = FALSE) {
+vstruct.detect = function(nodes, arcs, mb, data, alpha, extra.args = list(),
+    test, blacklist, max.sx = ncol(data), debug = FALSE) {
 
   vstruct.centered.on = function(x, mb, data, dsep.set) {
 
@@ -209,8 +204,8 @@ vstruct.detect = function(nodes, arcs, mb, data, alpha, B = NULL, test,
           cat("    > chosen d-separating set: '", sx, "'\n")
 
         a = allsubs.test(x = y, y = z, fixed = x, sx = sx, data = data,
-              test = test, B = B, alpha = alpha, max = min(max.sx, length(sx)),
-              debug = debug)
+              test = test, extra.args = extra.args, alpha = alpha,
+              max = min(max.sx, length(sx)), debug = debug)
 
         if (a["p.value"] <= alpha) {
 

R/backend-indep.R

@@ -1,14 +1,14 @@
 
 # do a single conditional independence test.
-ci.test = function(x, y, z, data, test, B, debug = FALSE) {
+ci.test = function(x, y, z, data, test, ..., debug = FALSE) {
 
   if (missing(x))
     stop("one or both of the variables to test are missing.")
 
   if (is.string(x)) {
 
-    ci.test.character(x = x, y = y, z = z, data = data, test = test, B = B,
-      debug = debug)
+    ci.test.character(x = x, y = y, z = z, data = data, test = test,
+      extra.args = list(...), debug = debug)
 
   }#THEN
   else if (is(x, c("matrix", "data.frame"))) {
@@ -22,7 +22,7 @@ ci.test = function(x, y, z, data, test, B, debug = FALSE) {
 
     nodes = names(x)
     ci.test.character(x = nodes[1], y = nodes[2], z = nodes[-(1:2)],
-      data = x, test = test, B = B, debug = debug)
+      data = x, test = test, extra.args = list(...), debug = debug)
 
   }#THEN
   else if (is.vector(x)) {
@@ -32,7 +32,7 @@ ci.test = function(x, y, z, data, test, B, debug = FALSE) {
 
     ci.test.vector(x = x, y = y, z = z, xlab = deparse(substitute(x)),
       ylab = deparse(substitute(y)), zlab = deparse(substitute(z)),
-      test = test, B = B, debug = debug)
+      test = test, extra.args = list(...), debug = debug)
 
   }#THEN
   else {
@@ -44,7 +44,7 @@ ci.test = function(x, y, z, data, test, B, debug = FALSE) {
 }#CI.TEST
 
 # do a single conditional independence test (nodes as character strings).
-ci.test.character = function(x, y, z, data, test, B, debug = FALSE) {
+ci.test.character = function(x, y, z, data, test, extra.args, debug = FALSE) {
 
   # the original data set is needed.
   data = check.data(data, allow.missing = TRUE, stop.if.all.missing = TRUE)
@@ -74,12 +74,13 @@ ci.test.character = function(x, y, z, data, test, B, debug = FALSE) {
   }#ELSE
   # check the test label.
   test = check.test(test, data = .data.frame.column(data, c(x, y, z)))
-  # check B (the number of permutation samples).
-  B = check.B(B, test)
+  # check the optional arguments to the test.
+  extra.args = check.test.args(test = test, extra.args = extra.args,
+                 data = .data.frame.column(data, c(x, y, z)))
 
   # create the htest object.
-  htest = indep.test(x = x, y = y, sx = z, data = data, test = test, B = B,
-            alpha = 1, learning = FALSE)
+  htest = indep.test(x = x, y = y, sx = z, data = data, test = test,
+            extra.args = extra.args, alpha = 1, learning = FALSE)
   htest$method = test.labels[test]
   htest$data.name = paste(x, "~", y, ifelse(length(z) > 0, "|", ""),
                       paste(z, collapse = " + "))
@@ -89,7 +90,8 @@ ci.test.character = function(x, y, z, data, test, B, debug = FALSE) {
 }#CI.TEST.CHARACTER
 
 # do a single conditional independence test (data vectors).
-ci.test.vector = function(x, y, z, xlab, ylab, zlab, test, B, debug = FALSE) {
+ci.test.vector = function(x, y, z, xlab, ylab, zlab, test, extra.args,
+    debug = FALSE) {
 
   # check debug.
   check.logical(debug)
@@ -141,12 +143,12 @@ ci.test.vector = function(x, y, z, xlab, ylab, zlab, test, B, debug = FALSE) {
   data = check.data(data, allow.missing = TRUE, stop.if.all.missing = TRUE)
   # check the test label.
   test = check.test(test, data = data)
-  # check B (the number of permutation samples).
-  B = check.B(B, test)
+  # check the optional arguments to the test.
+  extra.args = check.test.args(test = test, extra.args = extra.args, data = data)
 
   # create the htest object.
-  htest = indep.test(x = 1L, y = 2L, sx = sx, data = data, test = test, B = B,
-            alpha = 1, learning = FALSE)
+  htest = indep.test(x = 1L, y = 2L, sx = sx, data = data, test = test,
+            extra.args = extra.args, alpha = 1, learning = FALSE)
   htest$method = test.labels[test]
   htest$data.name = paste(xlab, "~", ylab,
         ifelse(length(z) > 0, paste("|", zlab), ""))

R/ci.test.R

@@ -1,6 +1,7 @@
 
 # reconstruct the equivalence class of a network.
-cpdag.backend = function(x, moral = FALSE, fix = FALSE, wlbl = TRUE, debug = FALSE) {
+cpdag.backend = function(x, moral = FALSE, fix = FALSE, wlbl = TRUE,
+    debug = FALSE) {
 
   nodes = names(x$nodes)
 

R/cpdag.R

@@ -434,6 +434,7 @@ weighting.sampling = function(fitted, event, evidence, n, batch, debug = FALSE)
     matching = r & !is.na(r)
 
     # compute the probabilities and use them as weigths.
+    attr(generated.data, "metadata") = collect.metadata(generated.data)
     w = weights(generated.data)
     cpe = cpe + sum(w[!is.na(r)])
     cpxe = cpxe + sum(w[matching])
@@ -460,7 +461,8 @@ weighting.sampling = function(fitted, event, evidence, n, batch, debug = FALSE)
 
 # generate random observations from conditional distributions with likelihood
 # weighting.
-weighting.distribution = function(fitted, nodes, evidence, n, batch, debug = FALSE) {
+weighting.distribution = function(fitted, nodes, evidence, n, batch,
+    debug = FALSE) {
 
   .Call(call_cpdist_lw,
         fitted = fitted,