R/initial_estimates.R

# Finds initial estimates for the parameters (R1, R2, C1) of the 1TM model.
# The estimates are inferred from the 3 estimated coefficients of a linear
# regression model by solving the non-linear equations
#
# Args:
#   flux    : Argument from one_tm()
#   data_df : Data frame containing the relevant data
#
# Returns: a numeric vector with names R1, R2 and C1
#
init_ests <- function(flux, data_df) {
  if (flux == "in") {
    in_init <- stats::lm(Qinp ~ x1 + x2 + Qinp1 - 1, data = data_df)
    sigma_hat <- sigma(in_init)
    b <- as.numeric(stats::coef(in_init))
  } else {
    out_init <- stats::lm(Qoutp ~ x1 + x3 + Qoutp1 - 1, data = data_df)
    sigma_hat <- sigma(out_init)
    b <- as.numeric(stats::coef(out_init))
  }
  tau <- data_df$tau[1]
  b1 <- b[1]
  b2 <- b[2]
  b3 <- b[3]
  r1 <- 2 * b2 / (b3 + 1)
  r2 <- r1 * b1 / (2 * r1 - b1 - 2 * b2)
  C1 <- tau * (r1 * r2 / b1 - (r1 + r2) / 2)
  R1 <- 1 / r1
  R2 <- 1 / r2
  if (flux == "in") {
    init <- c(R1 = R1, R2 = R2, C1 = C1)
  } else {
    init <- c(R1 = R2, R2 = R1, C1 = C1)
  }
  attr(init, "sigma_hat") <- sigma_hat
  return(init)
}
paulnorthrop/walls documentation built on May 15, 2019, 10:02 p.m.