R/2TM_regression_functions.R

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

# Functions using the parameterization (R1, R2, R3, C1, C2)

#' @export
two_tm_reg_fn <- function(x1, x2, x3, Qp1, R1, R2, R3, C1, C2, in_out, tau) {
  # Total number of Q values (Qin and Qout)
  n <- length(Qp1)
  m <- n / 2
  # Indices for inside and outside
#  ind_i <- which(in_out == "inside")
#  ind_o <- which(in_out == "outside")
  ind_i <- 1:m
  ind_o <- (m + 1): n
  # Quantities involved in the regression function
  r1 <- 1 / R1
  r2 <- 1 / R2
  r3 <- 1 / R3
  C1tau <- C1 / tau[1]
  C2tau <- C2 / tau[1]
  d1plus <- 2 * C1tau + r1 + r2
  d1minus <- 2 * C1tau - r1 - r2
  d2plus <- 2 * C2tau + r2 + r3
  d2minus <- 2 * C2tau - r2 - r3
  b1 <- d1minus / d1plus
  b2 <- d2minus / d2plus
  div <- 1 - r2 ^ 2 / (d1plus * d2plus)
  # Extract covariates for inside and for outside
  Qinp1 <- Qp1[ind_i]
  Qoutp1 <- Qp1[ind_o]
  x1i <- x1[ind_i]
  x2i <- x2[ind_i]
  x3i <- x3[ind_i]
  x1o <- x1[ind_o]
  x2o <- x2[ind_o]
  x3o <- x3[ind_o]
  # Qin regression function
  i1 <- 2 * r1 * r2 * (1 - r2 / d2plus) / d1plus
  i2 <- 2 * C1tau * r1 / d1plus
  i3 <- 2 * C2tau * r1 * r2 / (d1plus * d2plus)
  i4 <- (r2 ^ 2 / d2plus + d1minus) / d1plus
  i5 <- -r1 * r2 * (1 + d2minus / d2plus) / (r3 * d1plus)
  Qin_reg <- (i1 * x1i + i2 * x2i + i3 * x3i + i4 * Qinp1 + i5 * Qoutp1) / div
  # Qout regression function
  o1 <- 2 * r2 * r3 * (1 - r2 / d1plus) / d2plus
  o2 <- -2 * C1tau * r2 * r3 / (d1plus * d2plus)
  o3 <- -2 * C2tau * r3 / d2plus
  o4 <- -r2 * r3 * (1 + d1minus / d1plus) / (r1 * d2plus)
  o5 <- (r2 ^ 2 / d1plus + d2minus) / d2plus
  Qout_reg <- (o1 * x1o + o2 * x2o + o3 * x3o + o4 * Qinp1 + o5 * Qoutp1) / div
  return(c(Qin_reg, Qout_reg))
}

#' @export
two_tm_grad_fn <- function(x1, x2, x3, Qp1, R1, R2, R3, C1, C2, in_out, tau) {
  # R1, R2, R3, C1 and C2 are vectors when gnls() calls this, but not when
  # nls() calls.
  # Use R1[1], R2[1], R3[1], C1[1], C2[1]
  R1 <- R1[1]
  R2 <- R2[1]
  R3 <- R3[1]
  C1 <- C1[1]
  C2 <- C2[1]
  func <- function(x) {
    R1 <- x[1]
    R2 <- x[2]
    R3 <- x[3]
    C1 <- x[4]
    C2 <- x[5]
    val <- two_tm_reg_fn(x1, x2, x3, Qp1, R1, R2, R3, C1, C2, in_out, tau)
    return(val)
  }
  reg_fn <- two_tm_reg_fn(x1, x2, x3, Qp1, R1, R2, R3, C1, C2, in_out, tau)
  pars <- c(R1, R2, R3, C1, C2)
  for_jacobian <- list(func = func, x = pars)
  grad_mat <- do.call(numDeriv::jacobian, for_jacobian)
  dimnames(grad_mat) <- list(NULL, c("R1", "R2", "R3", "C1", "C2"))
  attr(reg_fn, "gradient") <- grad_mat
  return(reg_fn)
}

# Regression function using the parameterization
# (R1 + R2 + R3, R1 / (R1 + R2 + R3), R2 / (R1 + R2 + R3), C1, C2)

#' @export
R_two_tm_reg_fn <- function(x1, x2, x3, Qp1, sumR, fracR1, fracR2, C1, C2,
                            in_out, tau) {
  # Transform back to the parameterization (R1, R2, R3, C1, C2)
  R1 <- fracR1 * sumR
  R2 <- fracR2 * sumR
  R3 <- sumR - R1 - R2
  return(two_tm_reg_fn(x1, x2, x3, Qp1, R1, R2, R3, C1, C2, in_out, tau))
}

#' @export
R_two_tm_grad_fn <- function(x1, x2, x3, Qp1, sumR, fracR1, fracR2, C1, C2,
                             in_out, tau) {
  # sumR, fracR1, fracR2, C1 and C2 are vectors when gnls() calls this,
  # but not when nls() calls.
  # Use sumR[1], fracR1[1], fracR2[1], C1[1], C2[1]
  sumR <- sumR[1]
  fracR1 <- fracR1[1]
  fracR2 <- fracR2[1]
  C1 <- C1[1]
  C2 <- C2[1]
  func <- function(x) {
    sumR <- x[1]
    fracR1 <- x[2]
    fracR2 <- x[3]
    C1 <- x[4]
    C2 <- x[5]
    val <- R_two_tm_reg_fn(x1, x2, x3, Qp1, sumR, fracR1, fracR2, C1, C2,
                           in_out, tau)
    return(val)
  }
  reg_fn <- R_two_tm_reg_fn(x1, x2, x3, Qp1, sumR, fracR1, fracR2, C1, C2,
                            in_out, tau)
  pars <- c(sumR, fracR1, fracR2, C1, C2)
  for_jacobian <- list(func = func, x = pars)
  grad_mat <- do.call(numDeriv::jacobian, for_jacobian)
  dimnames(grad_mat) <- list(NULL, c("sumR", "fracR1", "fracR2", "C1", "C2"))
  attr(reg_fn, "gradient") <- grad_mat
  return(reg_fn)
}

# Regression function using the parameterization
# (U = 1 / (R1 + R2 + R3), R1 / (R1 + R2 + R3), R2 / (R1 + R2 + R3), C1, C2)

#' @export
U_two_tm_reg_fn <- function(x1, x2, x3, Qp1, Uvalue, fracR1, fracR2, C1, C2,
                            in_out, tau) {
  # Transform back to the parameterization (R1, R2, R3, C1, C2)
  sumR <- 1 / Uvalue
  R1 <- fracR1 * sumR
  R2 <- fracR2 * sumR
  R3 <- sumR - R1 - R2
  return(two_tm_reg_fn(x1, x2, x3, Qp1, R1, R2, R3, C1, C2, in_out, tau))
}

#' @export
U_two_tm_grad_fn <- function(x1, x2, x3, Qp1, Uvalue, fracR1, fracR2, C1, C2,
                             in_out, tau) {
  # sumR, fracR1, fracR2, C1 and C2 are vectors when gnls() calls this,
  # but not when nls() calls.
  # Use sumR[1], fracR1[1], fracR2[1], C1[1], C2[1]
  Uvalue <- Uvalue[1]
  fracR1 <- fracR1[1]
  fracR2 <- fracR2[1]
  C1 <- C1[1]
  C2 <- C2[1]
  func <- function(x) {
    Uvalue <- x[1]
    fracR1 <- x[2]
    fracR2 <- x[3]
    C1 <- x[4]
    C2 <- x[5]
    val <- U_two_tm_reg_fn(x1, x2, x3, Qp1, Uvalue, fracR1, fracR2, C1, C2,
                           in_out, tau)
    return(val)
  }
  reg_fn <- U_two_tm_reg_fn(x1, x2, x3, Qp1, Uvalue, fracR1, fracR2, C1, C2,
                            in_out, tau)
  pars <- c(Uvalue, fracR1, fracR2, C1, C2)
  for_jacobian <- list(func = func, x = pars)
  grad_mat <- do.call(numDeriv::jacobian, for_jacobian)
  dimnames(grad_mat) <- list(NULL, c("Uvalue", "fracR1", "fracR2", "C1", "C2"))
  attr(reg_fn, "gradient") <- grad_mat
  return(reg_fn)
}