Nothing
HiddenLayer <- function(number.of.nodes, prior = NULL,
expected.model.size = Inf) {
## Specify the structure of a hidden layer to be used in a feedforward neural
## network.
##
## TODO(steve): Include an option for a 'bias' term.
##
## Args:
## number.of.nodes: The number of output nodes in this hidden layer. The
## number of inputs is determined by the preceding layer.
## prior: An MvnPrior or SpikeSlabGlmPrior to use for the coefficients of
## each of the logistic regression models comprising this hidden layer.
## All models will use the same prior. The dimension of the prior must
## match the inputs from the previous layer.
##
## Returns:
## An object (list) encoding the necessary information for the underlying
## C++ code to build the desired neural network model.
check.scalar.integer(number.of.nodes)
stopifnot(number.of.nodes > 0)
stopifnot(is.null(prior) ||
inherits(prior, "MvnPrior") ||
inherits(prior, "SpikeSlabGlmPrior") ||
inherits(prior, "SpikeSlabGlmPriorDirect"))
ans <- list(number.of.nodes = number.of.nodes,
prior = prior,
expected.model.size = expected.model.size)
class(ans) <- c("HiddenLayerSpecification")
return(ans)
}
##===========================================================================
BayesNnet <- function(formula,
hidden.layers,
niter,
data,
subset,
prior = NULL,
expected.model.size = Inf,
drop.unused.levels = TRUE,
contrasts = NULL,
ping = niter / 10,
seed = NULL) {
## A Bayesian feed-forward neural network with logistic activation function
## and a Gaussian terminal layer.
##
## Args:
## formula: A model formula as one would pass to 'lm'.
## hidden.layers: A list of HiddenLayer objects defining the network
## structure.
## niter: The desired number of MCMC iterations.
## data: An optional data frame containing the variables used in 'formula'.
## subset: See 'lm'.
## prior: An object of class SpikeSlabPrior defining the prior distribution
## for the terminal layer. This includes the prior for the residual
## variance.
## drop.unused.levels: See 'lm'.
## contrasts: An optional list. See the 'contrasts.arg' of
## ‘model.matrix.default’.
## ping: The frequency with which to print status updates for the MCMC
## algorithm. Setting 'ping = 10' will print a status update message
## every 10 MCMC iterations.
## seed: The seed to use for the C++ random number generator.
##
## Returns:
## The MCMC draws for all the network coefficients. The return value also
## includes information needed by supporting methods (e.g. 'plot' and
## 'predict').
function.call <- match.call()
frame <- match.call(expand.dots = FALSE)
has.data <- !missing(data)
name.positions <- match(c("formula", "data", "subset", "na.action"),
names(frame), 0L)
frame <- frame[c(1L, name.positions)]
frame$drop.unused.levels <- drop.unused.levels
frame[[1L]] <- as.name("model.frame")
frame <- eval(frame, parent.frame())
model.terms <- attr(frame, "terms")
response <- model.response(frame, "numeric")
predictors <- model.matrix(model.terms, frame, contrasts)
## Check that each layer coheres with the preceding layer.
hidden.layers <- .CheckHiddenLayers(predictors, hidden.layers)
prior <- .EnsureTerminalLayerPrior(
response = response,
prior = prior,
hidden.layers = hidden.layers,
expected.model.size = expected.model.size)
check.positive.scalar(niter)
check.positive.scalar(ping)
if (!is.null(seed)) {
seed <- as.integer(seed)
}
ans <- .Call(analysis_common_r_do_feedforward,
predictors,
response,
hidden.layers,
prior,
as.integer(niter),
as.integer(ping),
seed)
ans$hidden.layer.specification <- hidden.layers
ans$niter <- niter
ans$contrasts <- attr(predictors, "contrasts")
ans$xlevels <- .getXlevels(model.terms, frame)
ans$call <- function.call
ans$terms <- model.terms
ans$response <- response
if (has.data) {
ans$training.data <- data
} else {
ans$training.data <- frame
}
dimnames(ans$hidden.layer.coefficients[[1]]) <-
list(NULL, colnames(predictors), NULL)
class(ans) <- c("BayesNnet")
return(ans)
}
##===========================================================================
plot.BayesNnet <- function(x, y = c("predicted", "residual", "structure",
"partial", "help"), ...) {
## Match the 'y' argument against the supplied default values, or the data frame
stopifnot(is.character(y))
which.function <- try(match.arg(y), silent = TRUE)
which.variable <- "all"
if (inherits(which.function, "try-error")) {
which.function <- "partial"
which.variable <- pmatch(y, names(x$training.data))
if (is.na(which.variable)) {
err <- paste0("The 'y' argument ", y, " must either match one of the plot types",
" or one of the predictor variable names.")
stop(err)
}
}
if (which.function == "predicted") {
PlotBayesNnetPredictions(x, ...)
} else if (which.function == "residual") {
PlotBayesNnetResiduals(x, ...)
} else if (which.function == "structure") {
PlotNetworkStructure(x, ...)
} else if (which.function == "help") {
help("plot.BayesNnet", package = "BoomSpikeSlab", help_type = "html")
} else if (which.function == "partial") {
if (which.variable == "all") {
varnames <- colnames(attributes(x$terms)$factors)
nvars <- length(varnames)
stopifnot(nvars >= 1)
nr <- max(1, floor(sqrt(nvars)))
nc <- ceiling(nvars / nr)
original.pars <- par(mfrow = c(nr, nc))
on.exit(par(original.pars))
for (i in 1:nvars) {
PartialDependencePlot(x, varnames[i], xlab = varnames[i], ...)
}
} else {
PartialDependencePlot(x, y, ...)
}
}
return(invisible(NULL))
}
##===========================================================================
SuggestBurn <- function(model) {
return(SuggestBurnLogLikelihood(-1 * model$residual.sd))
}
##===========================================================================
PlotBayesNnetPredictions <- function(model, burn = SuggestBurn(model), ...) {
pred <- predict(model, burn = burn)
predicted <- colMeans(pred)
actual <- model$response
plot(predicted, actual, ...)
abline(a = 0, b = 1)
}
##===========================================================================
PlotBayesNnetResiduals <- function(model, burn = SuggestBurn(model), ...) {
pred <- predict(model, burn = burn)
predicted <- colMeans(pred)
actual <- model$response
residual <- actual - predicted
plot(predicted, residual, ...)
abline(h = 0)
}
##===========================================================================
.HiddenLayerPriorMean <- function(prior) {
## Extract the mean of the prior distribution for hidden layer coefficients.
if (inherits(prior, "SpikeSlabPriorBase")) {
return(prior$mu)
} else if (inherits(prior, "MvnPrior")) {
return(prior$mean)
} else {
stop("Prior must be an MvnPrior or else inherit from SpikeSlabPriorBase")
}
}
##===========================================================================
.CheckHiddenLayers <- function(predictors, hidden.layers) {
## Check that the priors for each hidden layer are set and are of the
## appropriate dimension. Replace any NULL priors with default values.
## Return the updated list of hidden layers.
##
## Args:
## predictors: The matrix of predictors.
## hidden.layers: A list of HiddenLayerSpecification objects defining the
## hidden layers for the feed forward neural network model.
##
## Returns:
## hidden.layers, after checking that dimensions and priors are okay, and
## after replacing NULL priors with hopefully sensible defaults.
stopifnot(is.list(hidden.layers),
all(sapply(hidden.layers, inherits, "HiddenLayerSpecification")))
for (i in 1:length(hidden.layers)) {
if (is.null(hidden.layers[[i]]$prior)) {
expected.model.size <- hidden.layers[[i]]$expected.model.size
if (is.null(expected.model.size)) {
expected.model.size <- 10^10
}
check.positive.scalar(expected.model.size)
if (i == 1) {
## Logistic regression inputs are arbitrary.
prior <- LogitZellnerPrior(predictors,
prior.success.probability = .5,
expected.model.size = expected.model.size)
hidden.layers[[i]]$prior <- prior
} else {
## Logistic regression inputs are all between 0 and 1.
input.dimension <- hidden.layers[[i - 1]]$number.of.nodes
prior <- SpikeSlabGlmPriorDirect(
coefficient.mean = rep(0, input.dimension),
coefficient.precision = diag(rep(1, input.dimension)),
expected.model.size = expected.model.size)
hidden.layers[[i]]$prior <- prior
}
}
if (i == 1) {
stopifnot(length(.HiddenLayerPriorMean(hidden.layers[[i]]$prior))
== ncol(predictors))
} else {
stopifnot(length(.HiddenLayerPriorMean(hidden.layers[[i]]$prior))
== hidden.layers[[i - 1]]$number.of.nodes)
}
}
return(hidden.layers)
}
##===========================================================================
.EnsureTerminalLayerPrior <- function(response, hidden.layers, prior,
expected.model.size, ...) {
## Args:
## response: The vector of 'y' values from the regression.
## hidden.layers: A list of objects inheriting from
## HiddenLayerSpecification.
## prior: The prior distribution for the model in the terminal layer. This
## must be of type SpikeSlabPrior, SpikeSlabPriorDirect, or NULL. If NULL
## a default prior will be created.
##
## Returns:
## The checked prior distribution.
##
## Effects:
## Checks that the dimension of the prior matches the number of outputs in the
## final hidden layer.
dimension <- tail(hidden.layers, 1)[[1]]$number.of.nodes
if (is.null(prior)) {
precision <- diag(rep(1, dimension))
inclusion.probabilities <- rep(expected.model.size / dimension, dimension)
inclusion.probabilities[inclusion.probabilities > 1] <- 1
inclusion.probabilities[inclusion.probabilities < 0] <- 0
prior <- SpikeSlabPriorDirect(
coefficient.mean = rep(0, dimension),
coefficient.precision = precision,
prior.inclusion.probabilities = inclusion.probabilities,
sigma.guess = sd(response, na.rm = TRUE) / 2,
prior.df = 1)
}
stopifnot(inherits(prior, "SpikeSlabPrior")
|| inherits(prior, "SpikeSlabPriorDirect"))
return(prior)
}
##===========================================================================
predict.BayesNnet <- function(object, newdata = NULL, burn = 0, na.action = na.pass,
mean.only = FALSE, seed = NULL, ...) {
## Prediction method for BayesNnet.
## Args:
## object: object of class "BayesNnet" returned from the BayesNnet function.
## newdata: Either NULL, or else a data frame, matrix, or vector containing
## the predictors needed to make the prediction. If 'newdata' is 'NULL'
## then the predictors are taken from the training data used to create the
## model object. Note that 'object' does not store its training data, so
## the data objects used to fit the model must be present for the training
## data to be recreated. If 'newdata' is a data.frame it must contain
## variables with the same names as the data frame used to fit 'object'.
## If it is a matrix, it must have the same number of columns as
## object$beta. (An intercept term will be implicitly added if the number
## of columns is one too small.) If the dimension of object$beta is 1 or
## 2, then newdata can be a vector.
## burn: The number of MCMC iterations in 'object' that should be discarded.
## If burn <= 0 then all iterations are kept.
## na.action: what to do about NA's.
## mean.only: Logical. If TRUE then return the posterior mean of the
## predictive distribution. If FALSE then return the entire distribution.
## seed: Seed for the C++ random number generator.
## ...: extra aguments ultimately passed to model.matrix (in the event that
## newdata is a data frame)
## Returns:
## A matrix of predictions, with each row corresponding to a row in newdata,
## and each column to an MCMC iteration.
if (is.null(newdata)) {
predictor.matrix <- model.matrix(object, data = object$training.data)
} else {
predictor.matrix <- model.matrix(object, data = newdata)
}
stopifnot(is.matrix(predictor.matrix),
nrow(predictor.matrix) > 0)
check.nonnegative.scalar(burn)
check.scalar.boolean(mean.only)
if (!is.null(seed)) {
check.scalar.integer(seed)
}
ans <- .Call(analysis_common_r_feedforward_prediction,
object,
predictor.matrix,
as.integer(burn),
as.logical(mean.only),
seed)
class(ans) <- "BayesNnetPrediction"
return(ans)
}
##===========================================================================
PlotNetworkStructure <- function(model, ...) {
## Plot the nodes and edges of the neural network. Larger coefficients are
## thicker lines.
##
## Args:
## model: A model fit by BayesNnet.
## ...: Extra arguments passed to plot.igraph.
##
## NOTE: This function depends on having igraph installed. Igraph is a big
## package with lots of dependencies, so it is optional. If igraph is not
## installed then calling this function prints an error message.
ok <- requireNamespace("igraph", quietly = TRUE)
if (!ok) {
warning("Plotting network structure requires the 'igraph' package. ",
"Please install using 'install.packages'.")
return(NULL)
}
input.names <- dimnames(model$hidden.layer.coefficients[[1]])[[2]]
input.dimension <- length(input.names)
input.nodes <- data.frame(id = input.names,
layer = rep(0, length(input.names)),
position.in.layer = 1:length(input.names))
number.of.hidden.layers <- length(model$hidden.layer.specification)
hidden.node.counts <- sapply(model$hidden.layer.specification,
function(x) x$number.of.nodes)
## Edge weights are the absolute values of the coefficients in each layer,
## normalized by that layer so that input and output nodes don't dominate
## because of scaling issues.
layer <- rep(1:number.of.hidden.layers, times = hidden.node.counts)
position.in.layer <- c(sapply(hidden.node.counts, function(x) 1:x))
hidden.nodes <- data.frame(
id = paste("H", layer, position.in.layer, sep = "."),
layer = rep(1:number.of.hidden.layers, times = hidden.node.counts),
position.in.layer = position.in.layer)
terminal.node <- data.frame(
id = "terminal",
layer = length(hidden.node.counts) + 1,
position.in.layer = 1)
nodes <- rbind(input.nodes, hidden.nodes, terminal.node)
##---------------------------------------------------------------------------
## Compute the edges.
first.hidden.layer <- hidden.nodes[hidden.nodes$layer == 1, , drop = FALSE]
weights <- as.numeric(t(colMeans(model$hidden.layer.coefficients[[1]])))
nc <- max(abs(weights))
if (nc > 0) weights <- weights / nc
input.layer.edges <- data.frame(
from = as.character(rep(input.names, each = hidden.node.counts[1])),
to = as.character(rep(first.hidden.layer$id, times = length(input.names))),
weight = weights
)
edges <- input.layer.edges
if (number.of.hidden.layers >= 2) {
for (layer in 2:number.of.hidden.layers) {
current.layer <- nodes[nodes$layer == layer, , drop = FALSE]
previous.layer <- nodes[nodes$layer == layer - 1, , drop = FALSE]
weights <- as.numeric(t(colMeans(model$hidden.layer.coefficients[[layer]])))
nc <- max(abs(weights))
if (nc > 0) weights <- weights / nc
edges <- rbind(edges, data.frame(
from = as.character(rep(previous.layer$id, each = hidden.node.counts[layer])),
to = as.character(rep(current.layer$id, times = nrow(previous.layer))),
weight = weights
))
}
}
final.hidden.layer <-
nodes[nodes$layer == number.of.hidden.layers, , drop = FALSE]
weights <- colMeans(model$terminal.layer.coefficients)
nc <- max(abs(weights))
if (nc > 0) weights <- weights / nc
edges <- rbind(edges, data.frame(
from = as.character(final.hidden.layer$id),
to = as.character(terminal.node$id),
weight = weights
))
##---------------------------------------------------------------------------
## Compute the layout for the plot.
max.nodes <- max(nrow(input.nodes), hidden.node.counts)
initial.layer.node.offset <- (max.nodes - nrow(input.nodes)) / 2
hidden.layer.node.offsets <- (max.nodes - hidden.node.counts) / 2
terminal.layer.offset <- (max.nodes - 1) / 2
input.layer.layout <- cbind("layer" = 0,
"position.in.layer" = (1:nrow(input.nodes)) + initial.layer.node.offset)
hidden.layer.layout <- cbind(hidden.nodes[, c("layer", "position.in.layer")])
hidden.layer.layout[,2] <- hidden.layer.layout[, 2] +
hidden.layer.node.offsets[hidden.layer.layout[, 1]]
terminal.layer.layout <- cbind("layer" = 1 + number.of.hidden.layers,
"position.in.layer" = 1 + terminal.layer.offset)
graph.layout <- rbind(input.layer.layout, hidden.layer.layout,
terminal.layer.layout)
##---------------------------------------------------------------------------
## Do the plotting.
graph <- igraph::graph_from_data_frame(edges, vertices = NULL)
plot(graph, layout = as.matrix(graph.layout),
edge.color = edges$weight > 0,
edge.width = 5 * abs(edges$weight),
edge.arrow.size = 0,
...)
return(invisible(list(nodes = nodes, input.layer.edges = input.layer.edges)))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.