get_posterior_cov: Compute posterior covariance matrix block

In gschnabel/nucdataBaynet: Bayesian networks for nuclear data evaluation

Compute posterior covariance matrix block

Description

\loadmathjax

Compute the covariance matrix between selected elements according to the posterior distribution.

Usage

get_posterior_cov(map, zpost, U, obs, row_idcs, col_idcs, ret.dep = FALSE)

Arguments

map	Mapping object. Usually a compound map, see create_compound_map
zpost	Vector of posterior 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.
row_idcs	Indices associated with the variables of interest determining the number of rows of the posterior covariance matrix block
col_idcs	Indices associated with the variables of interest determining the number of columns of the posterior covariance matrix block
ret.dep	If TRUE, return the covariance matrix block for the dependent variables, otherwise of the independent variables. Default is FALSE.

Value

Return the posterior covariance matrix block \mjseqn\Sigma_I,J with the indices \mjseqnI specified in row_idcs and and the indices in \mjseqnJ specified in col_idcs.

Note

The posterior covariance matrix block is an approximation to the true posterior covariance matrix. It can be expected to be good if the non-linear relationships in the Bayesian network defined by map are in good approximation linear in the domain associated with significant values of the posterior density function.

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.