LMalgo: Maximum a posterior estimate via LM algorithm

In gschnabel/nucdataBaynet: Bayesian networks for nuclear data evaluation

Maximum a posterior estimate via LM algorithm

Description

Finds the values of the independent variables that maximize the value of the posterior probability density function using a customized Levenberg-Marquardt algorithm.

Usage

LMalgo(
  map,
  zprior,
  U,
  obs,
  zref = zprior,
  print.info = FALSE,
  adjust_idcs = NULL,
  ret.invcov = FALSE,
  control = list()
)

Arguments

map	Mapping object. Usually a compound map, see create_compound_map.
zprior	Vector of prior estimates of the independent variables (i.e., associated with nodes without parent nodes)
U	Prior covariance matrix of the independent variables
obs	Vector with observed values of dependent nodes. Must be of same size as zprior. An NA value in this vector means that the corresponding variable was not observed.
zref	Values of the independent variable that should be used as the initial reference point to calculate the Taylor approximation (using the jacobian function of map)
print.info	Display output regarding the progress of the optimization procedure
adjust_idcs	Indices of variables that should be adjusted, all other independent variables are assumed to be fixed. The default value NULL means that all independent variables are adjusted.
ret.invcov	If TRUE, also return the inverse posterior covariance matrix
control	A list with control parameters to tweak the convergence criteria, see Details.

Details

The convergence criteria of the Levenberg-Marquardt algorithm can be tweaked by several parameters in the control list passed as argument. The following fields are available:

reltol	If the relative improvement of posterior pdf value falls below this value for a number of subsequent iterations given by reltol_steps, terminate optimization. (Default 1e-6)
reltol_steps	Number of subsequent iteratios with a relative improvement of less than reltol necessary to terminate optimization procedure. (Default 3)
reltol2	If relative improvement falls below this value, immediately terminate the optimization procedure. (Default 1e-12)
mincount	Minimum number of iterations that should be performed regardless of relative improvement. (Default 10)
maxcount	Maximum number of iterations that should be performed regardless of relative improvement. (Default 100)
maxreject	Number of subsequent rejections of proposal vectors before the LM algorithm gives up. (Default 10)
tau	Parameter that influences the initial choice of the damping term. The larger this parameter, the larger the initial damping applied to the inverse posterior covariance matrix. (Default 1e-10)

Value

Return a list with the following fields:

zpost	Variable assignment corresponding to (potentially local) maximum of posterior density function.
init_val	Initial value of objective function (proportional to posterior pdf value)
final_val	Final value of objective function (proportional to posterior pdf value)
numiter	Number of iterations

Examples

library(Matrix)
params <- list(
  mapname = "mymap",
  maptype = "linearinterpol_map",
  src_idx = 1:10,
  tar_idx = 11:15,
  src_x = 1:10,
  tar_x = 3:7
)
mymap <- create_linearinterpol_map()
mymap$setup(params)

U <- Diagonal(n=15, x=c(rep(1e3, 10), rep(1, 5)))
zprior <- rep(0, 15)
zref <- rep(0, 15)
obs <- c(rep(NA,10), 5:9)

# glsalgo only works for linear relationships
# LMalgo can also deal with non-linear relationships
zpost <- glsalgo(mymap, zprior, U, obs)
optres <- LMalgo(mymap, zprior, U, obs)
zpost2 <- optres$zpost

# posterior estimates of values on computational grid
zpost[1:10]
# posterior estimates of error variables associated with observations
zpost[11:15]
# get posterior covariance block of independent variables
get_posterior_cov(mymap, zpost, U, obs, 1:5, 5:10)
# draw samples of independent variables from posterior distribution
get_posterior_sample(mymap, zpost, U, obs, 10)

gschnabel/nucdataBaynet documentation built on Feb. 3, 2023, 4:13 a.m.