R/averaged.fitted.R

In bnlearn: Bayesian Network Structure Learning, Parameter Learning and Inference

# average discrete nodes whose parameters are organised in a conditional
# probability table.
average.fitted.dnode = function(node, fitted, weights) {

  # create a conditional probability table filled with zeroes...
  cpt = fitted[[1]][[node]]$prob
  cpt[] = 0
  # ... add up the conditional probability tables from the networks, with
  # the respective weights...
  for (i in seq_along(fitted))
    cpt = cpt + fitted[[i]][[node]]$prob * weights[i]
  # ... and normalize the result so that columns sum up to 1 again.
  cpt = cpt / sum(weights)

  return(cpt)

}#AVERAGE.FITTED.DNODE

# average continuous nodes whose parameters are organised in vectors or matrices
# of regression coefficients and standard errors.
average.fitted.cgnode = function(node, fitted, weights) {

  # allocate the vector of the regression coefficients and the standard
  # error, both zeroed...
  coefs = fitted[[1]][[node]]$coefficients
  coefs[] = 0
  sd = fitted[[1]][[node]]$sd
  sd[] = 0
  # ... add up both of them, with the respective weights...
  for (i in seq_along(fitted)) {

    coefs = coefs + fitted[[i]][[node]]$coefficients * weights[i]
    sd = sd + fitted[[i]][[node]]$sd * weights[i]

  }#FOR
  # ... and normalize the result.
  coefs = coefs / sum(weights)
  sd = sd / sum(weights)

  return(list(coef = coefs, sd = sd))

}#AVERAGE.FITTED.CGNODE

# average of zero-inflated hyper-Poisson nodes.
average.fitted.zihpnode = function(node, fitted, weights) {

  # allocate the vectors of both sets of regression coefficients, zeroed...
  inflation = intensity = fitted[[1]][[node]]$inflation
  inflation[] = intensity[] = 0
  # ... and the dispersion parameter.
  dispersion = 0
  # add up all parameters, with the respective weights...
  for (i in seq_along(fitted)) {

    inflation = inflation + fitted[[i]][[node]]$inflation * weights[i]
    intensity = intensity + fitted[[i]][[node]]$intensity * weights[i]
    dispersion = dispersion + fitted[[i]][[node]]$dispersion * weights[i]

  }#FOR
  # ... and normalize the result.
  inflation = inflation / sum(weights)
  intensity = intensity / sum(weights)
  dispersion = dispersion / sum(weights)

  return(list(inflation = inflation, intensity = intensity,
              dispersion = dispersion))

}#AVERAGE.FITTED.ZIHPNODE

# average of zero-inflated negative binomial nodes.
average.fitted.zinbnode = function(node, fitted, weights) {

  # allocate the vectors of both sets of regression coefficients, zeroed...
  inflation = succprob = fitted[[1]][[node]]$inflation
  inflation[] = succprob[] = 0
  # ... and the number of failures.
  failures = 0
  # add up all parameters, with the respective weights...
  for (i in seq_along(fitted)) {

    inflation = inflation + fitted[[i]][[node]]$inflation * weights[i]
    succprob = succprob + fitted[[i]][[node]]$succprob * weights[i]
    failures = failures + fitted[[i]][[node]]$failures * weights[i]

  }#FOR
  # ... and normalize the result.
  inflation = inflation / sum(weights)
  succprob = succprob / sum(weights)
  failures = failures / sum(weights)

  return(list(inflation = inflation, succprob = succprob, failures = failures))

}#AVERAGE.FITTED.ZINBNODE

# average several bn.fit objects with the same structure.
average.fitted = function(fitted, weights) {

  # all the networks have the same structure and parameter sets, so we can
  # allocate the return value by copying one of them.
  averaged = fitted[[1]]
  cl = class(averaged)
  class(averaged) = "list"

  for (node in names(averaged)) {

    if (is(averaged[[node]], c("bn.fit.dnode", "bn.fit.onode"))) {

      averaged[[node]]$prob =
        average.fitted.dnode(node = node, fitted = fitted, weights = weights)

    }#THEN
    else if (is(averaged[[node]], c("bn.fit.gnode", "bn.fit.cgnode"))) {

      averaged[[node]][c("coefficients", "sd")] =
        average.fitted.cgnode(node = node, fitted = fitted, weights = weights)

      # in addition to averaging the parameters, remove the fitted values and
      # the residuals if present.
      averaged[[node]]$fitted.values = as.numeric(NA)
      averaged[[node]]$residuals = as.numeric(NA)
      # same with the configurations of the discrete parents in a conditional
      # Gaussian node.
      if ("configs" %in% names(averaged[[node]]))
        averaged[[node]]$configs =
          factor(NA, levels = levels(averaged[[node]]$configs))

    }#THEN
    else if (is(averaged[[node]], "bn.fit.zihpnode")) {

      averaged[[node]][c("inflation", "intensity", "dispersion")] =
        average.fitted.zihpnode(node = node, fitted = fitted, weights = weights)

      # in addition to averaging the parameters, remove the fitted values and
      # the residuals if present.
      averaged[[node]]$fitted.values = as.numeric(NA)
      averaged[[node]]$residuals = as.numeric(NA)

    }#THEN
    else if (is(averaged[[node]], "bn.fit.zinbnode")) {

      averaged[[node]][c("inflation", "succprob", "failures")] =
        average.fitted.zinbnode(node = node, fitted = fitted, weights = weights)

      # in addition to averaging the parameters, remove the fitted values and
      # the residuals if present.
      averaged[[node]]$fitted.values = as.numeric(NA)
      averaged[[node]]$residuals = as.numeric(NA)

    }#THEN

  }#FOR

  class(averaged) = cl

  # preserve attributes for classifiers when all fitted networks have the same
  # type and training variable.
  if (all(sapply(fitted, is, available.classifiers))) {

    classifier.type =
      sapply(fitted, function(x) intersect(class(x), available.classifiers))
    classifier.training = sapply(fitted, attr, "training")

    if ((length(unique(classifier.type)) == 1) &&
        (length(unique(classifier.training)) == 1)) {

      class(averaged) = union(unique(classifier.type), class(averaged))
      attr(averaged, "training") = unique(classifier.training)

    }#THEN

  }#THEN

  return(averaged)

}#AVERAGE.FITTED