NeuralEstimators: Likelihood-Free Parameter Estimation using Neural Networks

Description

An 'R' interface to the 'Julia' package 'NeuralEstimators.jl'. The package facilitates the user-friendly development of neural Bayes estimators, which are neural networks that map data to a point summary of the posterior distribution (Sainsbury-Dale et al., 2024, doi:10.1080/00031305.2023.2249522). These estimators are likelihood-free and amortised, in the sense that, once the neural networks are trained on simulated data, inference from observed data can be made in a fraction of the time required by conventional approaches. The package also supports amortised Bayesian or frequentist inference using neural networks that approximate the posterior or likelihood-to-evidence ratio (Zammit-Mangion et al., 2025, Sec. 3.2, 5.2, doi:10.48550/arXiv.2404.12484). The package accommodates any model for which simulation is feasible by allowing users to define models implicitly through simulated data.

Author(s)

Maintainer: Matthew Sainsbury-Dale msainsburydale@gmail.com

See Also

Useful links:

• https://github.com/msainsburydale/NeuralEstimators

• https://msainsburydale.github.io/NeuralEstimators.jl/dev/

assess a neural estimator

Description

assess a neural estimator

Usage

assess(
  estimators,
  parameters,
  Z,
  estimator_names = NULL,
  parameter_names = NULL,
  use_gpu = TRUE,
  verbose = TRUE
)

Arguments

Value

a list of two data frames: runtimes contains the total time taken for each estimator, while df is a long-form data frame with columns:

• "estimator"; the name of the estimator

• "parameter"; the name of the parameter

• "truth"; the true value of the parameter

• "estimate"; the estimated value of the parameter

• "m"; the sample size (number of iid replicates)

• "k"; the index of the parameter vector in the test set

• "j"; the index of the data set

See Also

risk(), rmse(), bias(), plotestimates(), and plotdistribution() for computing various empirical diagnostics and visualisations from an object returned by assess()

computes a Monte Carlo approximation of an estimator's bias

Description

computes a Monte Carlo approximation of an estimator's bias

Usage

bias(assessment, ...)

Arguments

Value

a dataframe giving the estimated bias

See Also

assess(), risk(), rmse()

bootstrap

Description

Compute bootstrap estimates from a neural point estimator

Usage

bootstrap(estimator, Z, B = 400, blocks = NULL, use_gpu = TRUE)

Arguments

Value

d × B matrix, where d is the dimension of the parameter vector

encodedata

Description

For data Z with missing (NA) entries, computes an augmented data set (U, W) where W encodes the missingness pattern as an indicator vector and U is the original data Z with missing entries replaced by a fixed constant c.

Usage

encodedata(Z, c = 0)

Arguments

Value

Augmented data set (U, W). If Z is provided as a list, the return type will be a JuliaProxy object; these objects can be indexed in the usual manner using [[, or converted to an R object using juliaGet() (note however that juliaGet() can be slow for large data sets).

See Also

the Julia version of encodedata()

Examples

## Not run: 
library("NeuralEstimators")
Z <- matrix(c(1, 2, NA, NA, 5, 6, 7, NA, 9), nrow = 3)
encodedata(Z)
encodedata(list(Z, Z))
## End(Not run)

estimate

Description

Apply a neural estimator to data

Usage

estimate(estimator, Z, X = NULL, batchsize = 32, use_gpu = TRUE)

Arguments

Value

a matrix of outputs resulting from applying estimator to Z (and possibly X)

See Also

sampleposterior() for making inference with neural posterior or likelihood-to-evidence-ratio estimators

Load a saved state of a neural estimator

Description

Load a saved state of a neural estimator (e.g., optimised neural-network parameters). Useful for amortised inference, whereby a neural network is trained once and then used repeatedly to make inference with new data sets.

Usage

loadstate(estimator, filename)

Arguments

Value

estimator updated with the saved state

Plot the empirical sampling distribution of an estimator.

Description

Plot the empirical sampling distribution of an estimator.

Usage

plotdistribution(
  df,
  type = c("box", "density", "scatter"),
  parameter_labels = NULL,
  estimator_labels = ggplot2::waiver(),
  truth_colour = "black",
  truth_size = 8,
  truth_line_size = NULL,
  pairs = FALSE,
  upper_triangle_plots = NULL,
  legend = TRUE,
  return_list = FALSE,
  flip = FALSE
)

Arguments

Value

a list of 'ggplot' objects or, if pairs = TRUE, a single 'ggplot'.

Examples

## Not run: 
# In the following, we have two estimators and, for each parameter, 50 estimates
# from each estimator.

estimators <- c("Estimator 1", "Estimator 2")
estimator_labels <- c("Estimator 1" = expression(hat(theta)[1]("·")),
                      "Estimator 2" = expression(hat(theta)[2]("·")))

# Single parameter:
df <- data.frame(
  estimator = estimators, truth = 0, parameter = "mu",
  estimate  = rnorm(2*50),
  replicate = rep(1:50, each = 2)
)
parameter_labels <- c("mu" = expression(mu))
plotdistribution(df)
plotdistribution(df, type = "density")
plotdistribution(df, parameter_labels = parameter_labels, estimator_labels = estimator_labels)

# Two parameters:
df <- rbind(df, data.frame(
  estimator = estimators, truth = 1, parameter = "sigma",
  estimate  = rgamma(2*50, shape = 1, rate = 1),
  replicate = rep(1:50, each = 2)
))
parameter_labels <- c(parameter_labels, "sigma" = expression(sigma))
plotdistribution(df, parameter_labels = parameter_labels)
plotdistribution(df, parameter_labels = parameter_labels, type = "density")
plotdistribution(df, parameter_labels = parameter_labels, type = "scatter")

# Three parameters:
df <- rbind(df, data.frame(
  estimator = estimators, truth = 0.25, parameter = "alpha",
  estimate  = 0.5 * runif(2*50),
  replicate = rep(1:50, each = 2)
))
parameter_labels <- c(parameter_labels, "alpha" = expression(alpha))
plotdistribution(df, parameter_labels = parameter_labels)
plotdistribution(df, parameter_labels = parameter_labels, type = "density")
plotdistribution(df, parameter_labels = parameter_labels, type = "scatter")
plotdistribution(df, parameter_labels = parameter_labels, type = "scatter", pairs = TRUE)

# Pairs plot with user-specified plots in the upper triangle:
upper_triangle_plots <- lapply(1:3, function(i) {
  x = rnorm(10)
  y = rnorm(10)
  shape = sample(c("Class 1", "Class 2"), 10, replace = TRUE)
  ggplot() +
    geom_point(aes(x = x, y = y, shape = shape)) + 
    labs(shape = "") +
    theme_bw()
})
plotdistribution(
    df, 
    parameter_labels = parameter_labels, estimator_labels = estimator_labels,
    type = "scatter", pairs = TRUE, upper_triangle_plots = upper_triangle_plots
    )
## End(Not run)

Plot estimates vs. true values.

Description

Plot estimates vs. true values.

Usage

plotestimates(
  df,
  estimator_labels = ggplot2::waiver(),
  parameter_labels = NULL
)

Arguments

Value

a 'ggplot' of the estimates for each parameter against the true value.

Examples

## Not run: 
K <- 50
df <- data.frame(
  estimator = c("Estimator 1", "Estimator 2"), 
  parameter = rep(c("mu", "sigma"), each = K),
  truth = 1:(2*K), 
  estimate = 1:(2*K) + rnorm(4*K)
)
estimator_labels <- c("Estimator 1" = expression(hat(theta)[1]("·")),
                      "Estimator 2" = expression(hat(theta)[2]("·")))
parameter_labels <- c("mu" = expression(mu), "sigma" = expression(sigma))

plotestimates(df,  parameter_labels = parameter_labels, estimator_labels)
## End(Not run)

posteriormode

Description

Computes the (approximate) posterior mode (maximum a posteriori estimate) given data Z.

Usage

posteriormode(estimator, Z, ...)

Arguments

Value

a d × K matrix of posterior samples, where d is the dimension of the parameter vector and K is the number of data sets provided in Z

See Also

sampleposterior() for sampling from the approximate posterior distribution

computes a Monte Carlo approximation of an estimator's Bayes risk

Description

computes a Monte Carlo approximation of an estimator's Bayes risk

Usage

risk(
  assessment,
  loss = function(x, y) abs(x - y),
  average_over_parameters = FALSE,
  average_over_sample_sizes = TRUE
)

Arguments

Value

a dataframe giving an estimate of the Bayes risk

See Also

assess(), bias(), rmse()

computes a Monte Carlo approximation of an estimator's root-mean-square error (RMSE)

Description

computes a Monte Carlo approximation of an estimator's root-mean-square error (RMSE)

Usage

rmse(assessment, ...)

Arguments

Value

a dataframe giving the estimated RMSE

See Also

assess(), bias(), risk()

sampleposterior

Description

Samples from the approximate posterior distribution given data Z.

Usage

sampleposterior(estimator, Z, N = 1000, ...)

Arguments

Value

a d × N matrix of posterior samples, where d is the dimension of the parameter vector. If Z is a list containing multiple data sets, a list of matrices will be returned

See Also

estimate() for making inference with neural Bayes estimators

save the state of a neural estimator

Description

save the state of a neural estimator

Usage

savestate(estimator, filename)

Arguments

Value

No return value, called for side effects

spatialgraph

Description

Constructs a graph object for use in a graph neural network (GNN).

Usage

spatialgraph(S, Z, isotropic = TRUE, stationary = TRUE, ...)

Arguments

Value

A GNNGraph (JuliaProxy object) or, if multiple data sets are provided, a vector of GNNGraph objects which can be indexed in the usual manner using [[ or converted to an R list using a combination of indexing and lapply.

Examples

## Not run: 
library("NeuralEstimators")

# Number of replicates
m <- 5

# Spatial locations fixed for all replicates
n <- 100 
S <- matrix(runif(n * 2), n, 2)
Z <- matrix(runif(n * m), n, m)
g <- spatialgraph(S, Z)

# Spatial locations varying between replicates
n <- sample(50:100, m, replace = TRUE)
S <- lapply(n, function(ni) matrix(runif(ni * 2), ni, 2))
Z <- lapply(n, function(ni) runif(ni))
g <- spatialgraph(S, Z)

# Multiple data sets: Spatial locations fixed for all replicates within a given data set
K <- 15 # number of data sets
n <- sample(50:100, K, replace = TRUE) # number of spatial locations can vary between data sets
S <- lapply(1:K, function(k) matrix(runif(n[k] * 2), n[k], 2))
Z <- lapply(1:K, function(k) matrix(runif(n[k] * m), n[k], m))
g <- spatialgraph(S, Z)

# Multiple data sets: Spatial locations varying between replicates within a given data set
S <- lapply(1:K, function(k) {
  lapply(1:m, function(i) {
  ni <- sample(50:100, 1)       # randomly generate the number of locations for each replicate
  matrix(runif(ni * 2), ni, 2)  # generate the spatial locations
  })
})
Z <- lapply(1:K, function(k) {
  lapply(1:m, function(i) {
    n <- nrow(S[[k]][[i]])
    runif(n)  
  })
})
g <- spatialgraph(S, Z)

## End(Not run)

tanhloss

Description

For k > 0, defines Julia code that defines the loss function,

L(\hat{\theta}, \theta) = \tanh\left(\frac{|\hat{\theta} - \theta|}{k}\right),

which approximates the 0-1 loss as k tends to zero.

The resulting string is intended to be used in the function train, but can also be converted to a callable function using juliaEval.

Usage

tanhloss(k)

Arguments

Value

String defining the tanh loss function in Julia code.

Train a neural estimator

Description

The function caters for different variants of "on-the-fly" simulation. Specifically, a sampler can be provided to continuously sample new parameter vectors from the prior, and a simulator can be provided to continuously simulate new data conditional on the parameters. If provided with specific sets of parameters (theta_train and theta_val) and/or data (Z_train and Z_val), they will be held fixed during training.

Note that using R functions to perform "on-the-fly" simulation requires the user to have installed the Julia package RCall.

Usage

train(
  estimator,
  sampler = NULL,
  simulator = NULL,
  theta_train = NULL,
  theta_val = NULL,
  Z_train = NULL,
  Z_val = NULL,
  m = NULL,
  M = NULL,
  K = 10000,
  xi = NULL,
  loss = "absolute-error",
  learning_rate = 1e-04,
  epochs = 100,
  batchsize = 32,
  savepath = NULL,
  stopping_epochs = 5,
  epochs_per_Z_refresh = 1,
  epochs_per_theta_refresh = 1,
  simulate_just_in_time = FALSE,
  use_gpu = TRUE,
  verbose = TRUE
)

Arguments

Value

a trained neural estimator or, if m is a vector, a list of trained neural estimators

See Also

assess() for assessing an estimator post training, and estimate()/sampleposterior() for making inference with observed data

Examples

## Not run: 
# Construct a neural Bayes estimator for replicated univariate Gaussian 
# data with unknown mean and standard deviation. 

# Load R and Julia packages
library("NeuralEstimators")
library("JuliaConnectoR")
juliaEval("using NeuralEstimators, Flux")

# Define the neural-network architecture
estimator <- juliaEval('
 n = 1    # dimension of each replicate
 d = 2    # number of parameters in the model
 w = 32   # width of each hidden layer
 psi = Chain(Dense(n, w, relu), Dense(w, w, relu))
 phi = Chain(Dense(w, w, relu), Dense(w, d))
 deepset = DeepSet(psi, phi)
 estimator = PointEstimator(deepset)
')

# Sampler from the prior
sampler <- function(K) {
  mu    <- rnorm(K)      # Gaussian prior for the mean
  sigma <- rgamma(K, 1)  # Gamma prior for the standard deviation
  theta <- matrix(c(mu, sigma), byrow = TRUE, ncol = K)
  return(theta)
}

# Data simulator
simulator <- function(theta_set, m) {
  apply(theta_set, 2, function(theta) {
    t(rnorm(m, theta[1], theta[2]))
  }, simplify = FALSE)
}

# Train using fixed parameter and data sets 
theta_train <- sampler(10000)
theta_val   <- sampler(2000)
m <- 30 # number of iid replicates
Z_train <- simulator(theta_train, m)
Z_val   <- simulator(theta_val, m)
estimator <- train(estimator, 
                   theta_train = theta_train, 
                   theta_val = theta_val, 
                   Z_train = Z_train, 
                   Z_val = Z_val)
                   
##### Simulation on-the-fly using R functions ####

juliaEval("using RCall") # requires the Julia package RCall
estimator <- train(estimator, sampler = sampler, simulator = simulator, m = m)

##### Simulation on-the-fly using Julia functions ####

# Defining the sampler and simulator in Julia can improve computational 
# efficiency by avoiding the overhead of communicating between R and Julia. 

juliaEval("using Distributions")

# Parameter sampler
sampler <- juliaEval("
      function sampler(K)
      	mu = rand(Normal(0, 1), K)
      	sigma = rand(Gamma(1), K)
      	theta = hcat(mu, sigma)'
      	return theta
      end")

# Data simulator
simulator <- juliaEval("
      function simulator(theta_matrix, m)
      	Z = [rand(Normal(theta[1], theta[2]), 1, m) for theta in eachcol(theta_matrix)]
      	return Z
      end")

# Train
estimator <- train(estimator, sampler = sampler, simulator = simulator, m = m)
## End(Not run)