R/simulation.R

In paulnorthrop/walls: Statistical Inference for the Thermophysical Properties of Walls from In-Situ Surveys

#' Simulate from the one thermal mass model (1TM)
#'
#' Simulates from the one thermal mass model (1TM)
#'
#' @param data A dataset in the format of either \code{\link{cwall_east}}
#'   or \code{\link{cwall_north}}.
#' @param tau A numeric scalar. Time interval between successive measurements,
#'   in hours.
#' @param sigma_in A numeric scalar.  Error standard deviation inside.
#' @param sigma_in A numeric scalar.  Error standard deviation outside.
#' @param n An integer scalar.
#' @param t0 A numeric scalar.  Initial temperature of the thermal mass C0.
#'   If \code{t0} is not supplied then the mean of the earliest
#' @examples
#' set.seed(23072018)
#' fit_n <- one_tm(data = cwall_north, flux = "both", hetero = FALSE)
#' params_n <- coef(fit_n)
#' params_n2 <- coef(fit_n2)
#' sigma(fit_n)
#'
#' # Estimates from Table 5.13 on page 173 of Gori (2017)
#' # Based on 3 days of data
#' pars <- c(R1 = 1.297, R2 = 0.069, C1 = 1.02e5 / 3600)
#' sim_n <- sim_1tm(data = cwall_north, params = pars, tc0 = 14.65,
#'                  start = "07/05/2015 14:00", end = "10/05/2015 13:55")
#' sim_n <- sim_1tm(data = cwall_north, params = pars, tc0 = 14.65,
#'                  start = "07/05/2015 14:00", end = "07/06/2015 13:55")
#' fit2 <- one_tm(data = cwall_north, flux = "both", sim_data = sim_n, hetero = FALSE)
#'
#' pjn <- sim_1tm(data = cwall_north, params = params_n, n = 1400)
#' pjn2 <- sim_1tm(data = cwall_north, params = params_n2, n = 1400)
#'
#' fit_n <- one_tm(data = cwall_north, flux = "both", hetero = FALSE)
#' params_n <- coef(fit_n)
#' sim_north <- sim_1tm(data = cwall_north, params = params_n, flux = "both")
#' sim_north <- sim_1tm(data = cwall_north, flux = "both")
#' fit_ns <- one_tm(data = cwall_north, flux = "both", sim_data = sim_north, hetero = FALSE)
#'
#' fit_e <- one_tm(data = cwall_east, flux = "both")
#' params_e <- coef(fit_e)
#' sim_east <- sim_1tm(data = cwall_east, params = params_e, flux = "both")
#' fit_es <- one_tm(data = cwall_east, flux = "out", sim_data = sim_east)
#' @seealso \code{\link{cwall_east}}, \code{\link{cwall_north}}.
sim_1tm <- function(data, start = NULL, end = NULL, params = NULL,
                    tau = 1 / 12, sigma_in = 1, sigma_out = 1, n = NULL,
                    tc0 = NULL) {
  if (!is.null(start) && !is.null(end)) {
    start_obs <- which(data[, "Time"] == start)
    end_obs <- which(data[, "Time"] == end)
    if (start_obs > end_obs) {
      stop("'start' must not be later than 'end'")
    }
    data <- data[start_obs:end_obs, ]
    n <- nrow(data)
  } else {
    if (is.null(n)) {
      n <- nrow(data)
    }
  }
  # Raw data
  # Check whether the format is like cwall_east or like cwall_north
  # Use internal and external surface temperatures, not air temperatures
  # This is partly because there tend to be fewer missing values in the
  # surface temperature data but also that the description of the prediction
  # of core temperatures on page 83 of Gori (2017) uses surface temperatures
  if ("T_in" %in% names(data)) {
    Qin <- data[1:n, "Q_in"]
    Qout <- data[1:n, "Q_out"]
    Ti <- data[1:n, "T_in"]
    Te <- data[1:n, "T_out"]
  } else {
    Qin <- data[1:n, "Q_in"]
    Qout <- data[1:n, "Q_out"]
    Ti <- data[1:n, "TS_int"]
    Te <- data[1:n, "TS_ext"]
  }
  #
  Ti_sum <- c(2 * Ti[1], Ti[2:n] + Ti[1:(n - 1)])
  Te_sum <- c(2 * Te[1], Te[2:n] + Te[1:(n - 1)])
  # Set default values of (R1, R2, C1) if params is not supplied
  if (is.null(params)) {
    params <- c(R1 = 1, R2 = 1, C1 = 100)
  }
  R1 <- params["R1"]
  R2 <- params["R2"]
  C1 <- params["C1"]
  #
  r1 <- 1 / R1
  r2 <- 1 / R2
  dminus <- 2 * C1 / tau - r1 - r2
  dplus <- 2 * C1 / tau + r1 + r2
  ap <- (Ti_sum / R1 + Te_sum / R2) / dplus
  b <- dminus / dplus
  #
  if (is.null(tc0)) {
    tc0 <- (Ti[1] + Te[1]) / 2
  }
  ############## To do: simulate errors in advance ... ##############
  sim_data <- data[1:n, ]
  sim_data[, "tc0"] <- rep(NA, n)
  sim_data[1, "tc0"] <- tc0
  sim_data[1, "Q_in"] <- (Ti[1] - tc0) / R1 +
    stats::rnorm(1, mean = 0, sd = sigma_in)
  sim_data[1, "Q_out"] <- (tc0 - Te[1]) / R2 +
    stats::rnorm(1, mean = 0, sd = sigma_out)
  for (i in 2:n) {
    tc0 <- ap[i] + b * tc0
    sim_data[i, "tc0"] <- tc0
    sim_data[i, "Q_in"] <- (Ti[i] - tc0) / R1 +
      stats::rnorm(1, mean = 0, sd = sigma_in)
    sim_data[i, "Q_out"] <- (tc0 - Te[i]) / R2 +
      stats::rnorm(1, mean = 0, sd = sigma_out)
  }
  #
  # NOTE: we negate the Qout values for consistency with the real datasets
  #
  sim_data[, "Q_out"] <- -sim_data[, "Q_out"]
  return(sim_data)
}