Nothing
# -------------------------------------------------------------------------------------------
#' Initialization of prior distributions
#'
#' @description Initializes the prior distributions for a model of class \code{TFPmodel} or
#' \code{NAWRUmodel}.
#'
#' @param model An object of class \code{TFPmodel} or \code{NAWRUmodel}.
#' @param MLE (Optional) A logical indicating whether the MLE estimates should be used for
#' the initialization. The default is \code{MLE = FALSE} if \code{MLEfit} is not provided
#' and vice versa.
#' @param MLEfit (Optional) An object of class \code{TFPfit} or \code{NAWRUfit} which is
#' used if \code{MLE = TRUE}.
#'
#' @return A list of three matrices with parameters for the prior distribution and box
#' constraints. Each list item refers to an equation, namely the \code{cycle}, \code{trend},
#' and second observation equation. Each list element is a \code{4 x n} matrix where \code{n}
#' denotes the number of parameters involved in the respective equation. The upper two
#' elements specify the distribution, the lower two parameters specify box constraints.
#' \code{NA} denotes no constraints. Autoregressive parameters are automatically restricted
#' to the stationary region unless box constraints are specified. The respective prior
#' distributions are defined through their mean and standard deviation. For instance,
#' \code{prior$cycle[, 1]} contains the mean, standard deviation, lower and upper bound for
#' the first variable, in that respective order.
#' @export
#'
#' @importFrom utils capture.output
initializePrior <- function(model, MLE = !is.null(MLEfit), MLEfit = NULL) {
# model class and attributes
class <- class(model)
trend <- attr(model, "trend")
cycle <- attr(model, "cycle")
attrib <- unlist(attributes(model), recursive = FALSE)
errorARMA <- attrib[grepl("errorARMA", names(attrib))][[1]]
# initialize list
prior <- list()
# load model data
namesExtract <- c("varName", "mean", "std", "LB", "UB")
prior <- .accessDfSystem(model = model)
prior <- lapply(prior, function(x) {
x[, namesExtract]
})
# rearrange data
prior <- lapply(prior, function(x) {
rownames(x) <- x[, "varName"]
x <- x[, c("mean", "std", "LB", "UB")]
x <- t(x)
return(x)
})
# drop trend variance
if (trend != "RW1") {
prior$trend <- prior$trend[, !grepl("tSigma", colnames(prior$trend)), drop = FALSE]
}
# obtain variance restrictions
varRestr <- .initializeVar(model = model, type = "hp", prior = TRUE, errorARMA = errorARMA)
trendNames <- colnames(prior$trend)[grepl("Sigma", colnames(prior$trend))]
# adjust variance priors
prior$cycle[c("mean", "std"), "cSigma"] <- mean(varRestr[, "cSigma"])
prior$trend[c("mean", "std"), trendNames] <- mean(varRestr[, 2:(dim(varRestr)[2] - 1)])
prior[[3]][c("mean", "std"), grepl("Sigma", colnames(prior[[3]]))] <- mean(varRestr[, "E2Sigma"])
# ----- priors from MLE
if (MLE) {
if (inherits(model, "TFPmodel")) {
if (!is.null(MLEfit) & !inherits(MLEfit, "TFPfit")) {
MLEfit <- NULL
message("The supplied object is not of class `TFPfit', starting MLE.")
}
if (is.null(MLEfit)) {
parRestr <- initializeRestr(model = model, type = "hp")
invisible(utils::capture.output(f <- suppressWarnings(fit(model = model, parRestr = parRestr))))
} else {
f <- MLEfit
}
pars <- f$parameters[sort(rownames(f$parameters)), 1]
names(pars) <- sort(rownames(f$parameters))
name <- rownames(f$parameters)
Tn <- length(f$tsl$tfpTrend)
# in case of NaNs for standard error; set standard error equal to estimate
f$parameters[is.na(f$parameters[, 2]), 2] <- f$parameters[is.na(f$parameters[, 2]), 1]
prior2 <- list()
prior2$cubs <- t(as.matrix(f$parameters[grepl("cu", name, fixed = TRUE), 1:2]))
prior2$trend <- t(as.matrix(f$parameters[grepl("td", name, fixed = TRUE), 1:2]))
prior2$cycle <- t(as.matrix(f$parameters[grepl("c", name, fixed = TRUE) & !grepl("cu", name, fixed = TRUE), 1:2]))
prior$cubs[1:2, colnames(prior2$cubs)] <- matrix(c(prior2$cubs[1, ], prior2$cubs[2, ] * sqrt(Tn)), 2, dim(prior2$cubs)[2], byrow = TRUE)
prior$trend[1:2, colnames(prior2$trend)] <- matrix(c(prior2$trend[1, ], prior2$trend[2, ] * sqrt(Tn)), 2, dim(prior2$trend)[2], byrow = TRUE)
prior$cycle[1:2, colnames(prior2$cycle)] <- matrix(c(prior2$cycle[1, ], prior2$cycle[2, ] * sqrt(Tn)), 2, dim(prior2$cycle)[2], byrow = TRUE)
prior$cubs[2, grepl("Sigma", colnames(prior$cubs))] <- prior$cubs[1, grepl("Sigma", colnames(prior$cubs))]
prior$trend[2, grepl("Sigma", colnames(prior$trend))] <- prior$trend[1, grepl("Sigma", colnames(prior$trend))]
prior$cycle[2, grepl("Sigma", colnames(prior$cycle))] <- prior$cycle[1, grepl("Sigma", colnames(prior$cycle))]
} else if (inherits(model, "NAWRUmodel")) {
if (!is.null(MLEfit) & !inherits(MLEfit, "NAWRUfit")) {
MLEfit <- NULL
message("The supplied object is not of class `NAWRUfit', starting MLE.")
}
if (is.null(MLEfit)) {
parRestr <- initializeRestr(model = model, type = "hp")
invisible(utils::capture.output(f <- suppressWarnings(.MLEfitNAWRU(model = model, parRestr = parRestr))))
} else {
f <- MLEfit
}
pars <- f$parameters[sort(rownames(f$parameters)), 1]
names(pars) <- sort(rownames(f$parameters))
name <- rownames(f$parameters)
Tn <- length(f$tsl$nawru)
# in case of NaNs for standard error; set standard error equal to estimate
f$parameters[is.na(f$parameters[, 2]), 2] <- f$parameters[is.na(f$parameters[, 2]), 1]
prior2 <- list()
prior2$pcInd <- t(as.matrix(f$parameters[grepl("pc", name, fixed = TRUE), 1:2]))
prior2$trend <- t(as.matrix(f$parameters[grepl("td", name, fixed = TRUE), 1:2]))
prior2$cycle <- t(as.matrix(f$parameters[grepl("c", name, fixed = TRUE) & !grepl("pc", name, fixed = TRUE), 1:2]))
prior$pcInd[1:2, colnames(prior2$pcInd)] <- matrix(c(prior2$pcInd[1, ], prior2$pcInd[2, ] * sqrt(Tn)), 2, dim(prior2$pcInd)[2], byrow = TRUE)
prior$trend[1:2, colnames(prior2$trend)] <- matrix(c(prior2$trend[1, ], prior2$trend[2, ] * sqrt(Tn)), 2, dim(prior2$trend)[2], byrow = TRUE)
prior$cycle[1:2, colnames(prior2$cycle)] <- matrix(c(prior2$cycle[1, ], prior2$cycle[2, ] * sqrt(Tn)), 2, dim(prior2$cycle)[2], byrow = TRUE)
prior$pcInd[2, grepl("Sigma", colnames(prior$pcInd))] <- prior$pcInd[1, grepl("Sigma", colnames(prior$pcInd))]
prior$trend[2, grepl("Sigma", colnames(prior$trend))] <- prior$trend[1, grepl("Sigma", colnames(prior$trend))]
prior$cycle[2, grepl("Sigma", colnames(prior$cycle))] <- prior$cycle[1, grepl("Sigma", colnames(prior$cycle))]
}
}
# ensure unimodal distribution for beta for RAR2 cycle (alpha, beta (shapes) > 1)
if (cycle == "RAR2") {
m1 <- (prior$cycle[1, 1] - prior$cycle[3, 1]) / (prior$cycle[4, 1] - prior$cycle[3, 1])
m2 <- prior$cycle[2, 1]
cond1 <- sqrt(m1^2 * (1 - m1) / (1 + m1)) * (prior$cycle[4, 1] - prior$cycle[3, 1])
cond2 <- sqrt((m1^3 - 2 * m1^2 + m1) / (2 + m1)) * (prior$cycle[4, 1] - prior$cycle[3, 1])
prior$cycle[2, 1] <- min(cond1, cond2, m2)
m1 <- (prior$cycle[1, 2] - prior$cycle[3, 2]) / (prior$cycle[4, 2] - prior$cycle[3, 2])
m2 <- prior$cycle[2, 2]
cond1 <- sqrt(m1^2 * (1 - m1) / (1 + m1)) * (prior$cycle[4, 2] - prior$cycle[3, 2])
cond2 <- sqrt((m1^3 - 2 * m1^2 + m1) / (2 + m1)) * (prior$cycle[4, 2] - prior$cycle[3, 2])
prior$cycle[2, 2] <- min(cond1, cond2, m2)
}
# ensure that standard deviation is not too small (invertible)
sd_min <- 1e-6 # such that var_min = 1e-12
prior[[1]][2, ] <- sapply(prior[[1]][2, ], max, sd_min)
prior[[2]][2, ] <- sapply(prior[[2]][2, ], max, sd_min)
prior[[3]][2, ] <- sapply(prior[[3]][2, ], max, sd_min)
# return
return(prior)
}
# -------------------------------------------------------------------------------------------
#' Converts the mean and standard deviation of a Gamma-distributed variable into the
#' parameters \eqn{s} and \eqn{\nu}.
#'
#' @param m mean parameter
#' @param std standard deviation
#'
#' @details The parameters \code{s} and \code{nu} are related to the regular shape
#' \eqn{\alpha} and rate \eqn{\beta} parametrization in the following way:
#' \eqn{\alpha = \nu / 2}
#' \eqn{\beta = s / 2}
#'
#' @return A vector with parameters \eqn{s} and \eqn{\nu}.
#' @keywords internal
.meanStd2GammasNu <- function(m, std) {
nu <- 2 * m^2 / std^2 + 4
s <- m * (nu - 2)
result <- c(s, nu)
return(result)
}
# -------------------------------------------------------------------------------------------
#' Transforms the prior distribution defined by mean and standard deviation to the
#' appropriate input parameters.
#'
#' @param prior A matrix containing the mean and standard deviation of the priors.
#' @param restr A matrix containing box constraints for the parameters.
#' @param namesInvGammaDistr A character vector specifying the names of all inverse Gamma
#' distributed parameters.
#' @param namesNormalDistr A character vector specifying the names of all normally
#' distributed parameters.
#' @param namesBetaDistr A character vector specifying the names of all Beta distributed
#' parameters.
#' @keywords internal
.priorMSd2Parameter <- function(prior, restr, namesInvGammaDistr, namesNormalDistr, namesBetaDistr) {
# transform mean and standard deviation to appropriate input parameters
for (names in namesInvGammaDistr) {
prior[, names] <- .meanStd2GammasNu(m = prior[1, names], std = prior[2, names])
}
for (names in namesBetaDistr) {
prior[, names] <- .meanStd2Beta(
m = prior[1, names], std = prior[2, names],
lb = restr[1, names], ub = restr[2, names]
)
}
prior[2, namesNormalDistr] <- prior[2, namesNormalDistr]^2
# deal with no box constraints
distrPar <- rbind(prior, restr)
distrPar[3, is.na(distrPar[3, ])] <- -Inf
distrPar[4, is.na(distrPar[4, ])] <- Inf
# return
return(distrPar)
}
# -------------------------------------------------------------------------------------------
#' Converts the mean and standard deviation of a (possibly scaled) Beta-distributed variable
#' into the two shape parameters \eqn{\alpha} and \eqn{\beta}.
#'
#' @param m mean parameter
#' @param std standard deviation
#' @param lb lower bound
#' @param ub upper bound
#'
#' @return A vector with shape parameters \eqn{\alpha} and \eqn{\beta}.
#' @keywords internal
.meanStd2Beta <- function(m, std, lb = 0, ub = 1) {
# unscaled mean and standard deviation
m <- (m - lb) / (ub - lb)
std <- std / (ub - lb)
# check bounds
if (m < 0 | m > 1) {
m <- (m < 0) * (0 + 1e-8) + (m > 1) * (1 - 1e-8)
}
if (std^2 > m * (1 - m)) {
std <- sqrt(m * (1 - m) - 1e-8)
}
alpha <- m^2 / std^2 * (1 - m) - m
beta <- alpha / m - alpha
result <- c(alpha, beta)
return(result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.