R/bme_estimate.R

In BMEmapping: Spatial Interpolation using Bayesian Maximum Entropy (BME)

# ============================================================================
# Wrapper to bme_mean and bme_mode functions
#
# Details:
# Compute the bme mean, variance and mode of an estimation location
#
# Inputs:
# -  x matrix of estimation locations
# -  ch matrix of hard data locations
# -  cs matrix of soft data locations
# -  zh vector of hard data
# -  a vector of lower bounds of soft data
# -  b vector of lower bounds of soft data
# -  model string name of covariance or variogram model
# -  nugget a non-negative value
# -  sill a non-negative value
# -  range a non-negative value
# -  nsmax number of soft data locations closer to the estimation location
# -  nhmax number of hard data locations closer to the estimation location
#
# Outputs:
# - A data frame of estimation locations with their corresponding bme mean,
#   variance and mode estimates
# ============================================================================
bme_estimate <- function(x, ch, cs, zh, a, b, model, nugget, sill, range,
                         nsmax = 5, nhmax = 5, n = 50,
                         zk_range = extended_range(zh, a, b)) {
  #x <- matrix(c(x), ncol = 2)
  nk <- nrow(x)

  # set up container for estimates: mean, variance, mode
  df <- matrix(NA, ncol = 3, nrow = nk)

  for (i in 1:nk) {
    d <- prob_zk(
      x = x[i, ], ch = ch, cs = cs, zh = zh, a = a, b = b,
      model = model, nugget = nugget, sill = sill, range = range,
      nsmax = nsmax, nhmax = nhmax, n = n, zk_range = zk_range
    )

    delta <- d[2, 1] - d[1, 1]

    # compute mean
    zk_mean <- sum(d[, 1] * d[, 2] * delta)

    # compute variance
    zk_var <- sum((d[, 1] - zk_mean)^2 * d[, 2] * delta)

    # compute mode
    zk_mode <- d[which.max(d[, 2]), 1]

    # gather estimates
    df[i, ] <- round(c(zk_mode, zk_mean, zk_var), 4)
  }

  return(df)
}