R/priors.R

In ashiklom/edr-da: ED2 RTM testing functions

#' Multivariate allometry prior (varying both base and exponent)
#'
#' @param traits Traits list
#' @export
rallom2 <- function() {
  purrr::map2(
    allom_mu,
    allom_Sigma,
    ~mvtnorm::rmvnorm(1, .x, .y)[1, ]
  ) %>%
    purrr::map(~c(exp(.[1]), .[2]))
}

#' @rdname rallom2
#' @export
dallom2 <- function(traits, log = TRUE) {
  b1Bl <- purrr::map_dbl(traits, allom_names[1]) %>% log()
  b2Bl <- purrr::map_dbl(traits, allom_names[2])
  pfts <- names(traits)
  purrr::pmap_dbl(
    list(b1Bl, b2Bl, pfts),
    ~mvtnorm::dmvnorm(c(..1, ..2), allom_mu[[..3]], allom_Sigma[[..3]], log = log)
  )
}

#' @rdname rallom2
#' @export
rwallom2 <- function() {
  purrr::map2(
    wallom_mu,
    wallom_Sigma,
    ~mvtnorm::rmvnorm(1, .x, .y)[1, ]
  ) %>%
    purrr::map(~c(exp(.[1]), .[2]))
}

#' @rdname rallom2
#' @export
dwallom2 <- function(traits, log = TRUE) {
  b1Bw <- purrr::map_dbl(traits, wallom_names[1]) %>% log()
  b2Bw <- purrr::map_dbl(traits, wallom_names[2])
  pfts <- names(traits)
  purrr::pmap_dbl(
    list(b1Bw, b2Bw, pfts),
    ~mvtnorm::dmvnorm(c(..1, ..2), wallom_mu[[..3]], wallom_Sigma[[..3]], log = log)
  )
}

#' Allometry mean vector
#'
#' @name allom_mu
#' @export allom_mu
NULL

#' Wood allometry mean vector
#'
#' @name wallom_mu
#' @export wallom_mu
NULL

#' Allometry variance-covariance matrix
#'
#' @name allom_Sigma
#' @export allom_Sigma
NULL

#' Wood allometry variance-covariance matrix
#'
#' @name wallom_Sigma
#' @export wallom_Sigma
NULL

#' Univariate allometry prior (varying just the base)
#'
#' @inheritParams rallom2
#' @export
rallom1 <- function() {
  out <- rnorm(npfts, b1Bl_means, b1Bl_sds) %>% exp()
  names(out) <- paste(pfts, allom_names[1], sep = ".")
  out
}

#' @rdname rallom1
#' @export
dallom1 <- function(traits, log = TRUE) {
  allom_vals <- purrr::map_dbl(traits, allom_names[1]) %>% log()
  ld <- dnorm(allom_vals, b1Bl_means, b1Bl_sds, log = log)
  ld
}

#' @rdname rallom1
#' @export
rwallom1 <- function() {
  out <- rnorm(npfts, b1Bw_means, b1Bw_sds) %>% exp()
  names(out) <- paste(pfts, wallom_names[1], sep = ".")
  out
}

#' @rdname rallom1
#' @export
dwallom1 <- function(traits, log = TRUE) {
  allom_vals <- purrr::map_dbl(traits, wallom_names[1]) %>% log()
  ld <- dnorm(allom_vals, b1Bw_means, b1Bw_sds, log = log)
  ld
}


#' @export
prior_clumping <- c(0, 1)

#' Priors on other EDR parameters
#'
#' @inheritParams rallom2
#' @export
rclumping <- function() {
  out <- runif(npfts, prior_clumping[1], prior_clumping[2])
  names(out) <- paste(pfts, "clumping_factor", sep = ".")
  out
}

#' @rdname rclumping
#' @export
dclumping <- function(traits, log = TRUE) {
  x <- purrr::map_dbl(traits, "clumping_factor")
  out <- dunif(x, prior_clumping[1], prior_clumping[2], log = log)
  names(out) <- pfts
  out
}

#' @export
prior_orient <- c(6, 4) * 3

#' @rdname rclumping
#' @export
rorient <- function() {
  ntry <- 1
  out <- 2 * rbeta(length(pfts), prior_orient[1], prior_orient[2]) - 1
  while (any(out > 0.6)) {
    ibad <- out > 0.6
    ntry <- ntry + 1
    if (ntry > 100) stop("Too many tries. Something is wrong with orient.")
    out[ibad] <- 2 * rbeta(sum(ibad), prior_orient[1], prior_orient[2]) - 1
  }
  names(out) <- paste(pfts, allom_names[1], sep = ".")
  out
}

#' @rdname rclumping
#' @export
dorient <- function(traits, log = TRUE) {
  x <- purrr::map_dbl(traits, "orient_factor")
  out <- dbeta((x + 1) / 2, prior_orient[1], prior_orient[2], log = log)
  names(out) <- pfts
  out
}

#' @export
prior_residual <- c(0.03, 1.5)

#' @rdname rclumping
#' @export
rresidual <- function(n = 1) {
  rgamma(n, prior_residual[1], prior_residual[2]) %>%
    setNames(paste0("residual", seq_len(n)))
}

#' @rdname rclumping
#' @export
dresidual <- function(params, log = TRUE) {
  if (!is.null(names(1:5))) {
    x <- params["residual"]
  } else {
    # HACK: Assume it's the last value
    x <- tail(params, 1)
  }
  dgamma(x, prior_residual[1], prior_residual[2], log = log) %>%
    setNames(paste0("residual", seq_len(n)))
}

#' @export
prior_residual2 <- list(
  intercept = 10,
  slope = 1
)

#' @rdname rclumping
#' @export
rresidual2 <- function(n = 1) {
  rint <- rexp(n, prior_residual2$intercept)
  names(rint) <- paste0("residual_intercept", seq_len(n))
  rslope <- rexp(n, prior_residual2$slope)
  names(rslope) <- paste0("residual_slope", seq_len(n))
  c(rint, rslope)
}

#' @rdname rclumping
#' @export
dresidual2 <- function(params, log = TRUE) {
  rint <- dnorm(
    params["residual_intercept"],
    prior_residual2$intercept[1],
    prior_residual2$intercept[2],
    log = log
  )
  rslope <- dnorm(
    params["residual_slope"],
    prior_residual2$slope[1],
    prior_residual2$slope[2],
    log = log
  )
  c("residual_intercept" = rint, "residual_slope" = rslope)
}

#' @rdname rclumping
#' @export
rprospect <- function() {
  out <- purrr::map(
    pfts,
    ~rmvnorm_positive(prospect_means[., ], prospect_covars[, , .])
  )
  names(out) <- pfts
  out
}

#' @rdname rclumping
#' @export
dprospect <- function(traits, log = TRUE) {
  out <- purrr::imap_dbl(
    traits,
    ~mvtnorm::dmvnorm(.x[prospect_names], prospect_means[.y, ], prospect_covars[, , .y], log = TRUE)
  )
  names(out) <- pfts
  out
}

#' Draw only positive values from multivariate prior
rmvnorm_positive <- function(mu, Sigma) {
  draw <- -1
  while(any(draw < 0)) {
    draw <- mvtnorm::rmvnorm(1, mu, Sigma)[1,]
  }
  draw
}