R/model_functions.R

In BayerSe/trexreg: Truncated Expectation Regression

#' @title Specification Function
#' @description TBA
#' @param x numeric vectors
#' @keywords internal
#' @name spec_fun
NULL

#' @rdname spec_fun
g <- function(x) 0

#' @rdname spec_fun
g_prime <- function(x) 0

#' @rdname spec_fun
g1 <- function(x, a, b) 2/(b-a) * x

#' @rdname spec_fun
g1_prime <- function(x, a, b) 2/(b-a) * 1

#' @rdname spec_fun
g2 <- function(x) x

#' @rdname spec_fun
g2_prime <- function(x) 1

#' @rdname spec_fun
phi2 <- function(x) -1/x

#' @rdname spec_fun
phi2_prime <- function(x) 1/x^2

#' @rdname spec_fun
phi <- function(x) log1p(exp(x))

#' @rdname spec_fun
phi_prime <- function(x) 1 / (1 + exp(-x))

#' @rdname spec_fun
phi_prime2 <- function(x)  exp(-x) / (1 + exp(-x))^2


#' @title Loss function
#' @description TBA
#' @inherit psi_e
#' @return TBA
#' @keywords internal
#' @export
rho_e <- function(y, qa, ea, uq) {
  idx <- (y <= qa) * 1
  v1 <- (idx - uq) * g(qa) - idx * g(y)
  v2 <- phi2_prime(ea) * (ea - qa + (qa - y) * idx / uq) - phi2(ea)
  v <- v1 + v2
  mean(v)
}

#' @title Loss function
#' @description TBA
#' @inherit psi_m
#' @return TBA
#' @keywords internal
#' @export
rho_m <- function(y, qa, qb, mab, lq, uq) {
  idx_lq <- (y <= qa) * 1
  idx_uq <- (y <= qb) * 1
  v1 <- (idx_lq - lq) * g1(qa, lq, uq) - idx_lq * g1(y, lq, uq)
  v2 <- (idx_uq - uq) * g2(qb) - idx_uq * g2(y)
  v3 <- phi_prime(mab) *
    (mab + 1/(uq-lq) * (idx_uq * (qb - y) - uq*qb - idx_lq * (qa - y) + lq*qa)) -
    phi(mab)

  mean(v1 + v2 + v3)
}

#' @title Regression model loss function
#' @description TBA
#' @param par TBA
#' @inherit lter.fit
#' @return TBA
#' @keywords internal
#' @export
uter.loss <- function(par, xm, xu, y, uq) {
  k1 <- ncol(xm)
  k2 <- ncol(xu)
  par_xm <- par[1:k1]
  par_xu <- par[(k1+1):(k1+k2)]
  ea <- ter.lin_mod(par = par_xm, x = xu)
  qa <- ter.lin_mod(par = par_xu, x = xm)
  v <- rho_e(y = y, qa = qa, ea = ea, uq = uq)
  mean(v)
}

#' @title Regression model loss function
#' @description TBA
#' @param par TBA
#' @inherit dter.fit
#' @return TBA
#' @keywords internal
#' @export
dter.loss <- function(par, xm, xl, xu, y, lq, uq) {
  k1 <- ncol(xm)
  k2 <- ncol(xl)
  k3 <- ncol(xu)

  par_m <- par[1:k1]
  par_qa <- par[(k1+1):(k1+k2)]
  par_qb <- par[(k1+k2+1):(k1+k2+k3)]

  mab <- ter.lin_mod(par = par_m, x = xm)
  qa <- ter.lin_mod(par = par_qa, x = xl)
  qb <- ter.lin_mod(par = par_qb, x = xu)

  v <- rho_m(y = y, qa = qa, qb = qb, mab = mab, lq = lq, uq = uq)
  mean(v)
}

#' @title Identification function
#' @description TBA
#' @param y TBA
#' @param qa TBA
#' @param gqa TBA
#' @param ea TBA
#' @param gea TBA
#' @param return_mean TBA
#' @inherit ster.fit
#' @return TBA
#' @keywords internal
#' @export
psi_e <- function(y, qb, gqb, gmb, mb, uq, return_mean=TRUE) {
  k <- ncol(gqb)
  idx <- (y <= qb) * 1
  v1 <- gqb * replicate(k, as.numeric((idx - uq) / uq * (-1/mb)))
  v2 <- gmb * replicate(k, as.numeric(1 / mb^2 * (mb - qb + (qb - y) * idx / uq)))
  v <- v1 + v2
  if (return_mean) {
    colMeans(v)
  } else {
    v
  }
}

#' @title Loss function
#' @description TBA
#' @param y TBA
#' @param qa TBA
#' @param qb TBA
#' @param gqa TBA
#' @param gqb TBA
#' @param mab TBA
#' @param qmab TBA
#' @param return_mean TBA
#' @inherit dter.fit
#' @return TBA
#' @keywords internal
#' @export
psi_m <- function(y, qa, qb, mab, gqa, gqb, gmab, lq, uq, return_mean=TRUE) {
  k <- ncol(gqa)
  idx_lq <- (y <= qa) * 1
  idx_uq <- (y <= qb) * 1

  v1 <- gqa * replicate(k, as.numeric((idx_lq - lq) * (g1_prime(qa, lq, uq) - phi_prime(mab) / (uq-lq))))
  v2 <- gqb * replicate(k, as.numeric((idx_uq - uq) * (g2_prime(qb) + phi_prime(mab) / (uq-lq))))
  v3 <- gmab * replicate(k, as.numeric(phi_prime2(mab) *
                                         (mab + 1/(uq-lq) * (idx_uq * (qb - y) - uq*qb - idx_lq  * (qa - y) + lq*qa))))
  v <- v1 + v2 + v3
  if (return_mean) {
    colMeans(v)
  } else {
    v
  }
}