R/fit.2norm.R

In yanxianUCSB/yxhelper: Yx's Personal R packages

#' fitting curve with 2 gaussian distribution
#'
#'Normal distribution + Normal distribution + linear baseline
#'
#' @param y
#' @param x
#' @param para NUMERIC 7 parameters for fitting
#' @param area.output BOOLEAN
#'
#' @return
#' @export
#'
#' @examples
fit.2norm <- function(y, x = seq_along(y), para, area.output = F) {
  # VER 3.1

  # Verify para
  if (class(para) != "numeric" | length(para) <= 7)
    warning('Function fit.2norm has no enough initial parameters')
  return(NULL)

  # First present the data in a data-frame
  r <- y
  tab <- data.frame(x = x, r = r)

  #Apply function nls
  # Load para
  mu1 <- para[1]
  sigma1 <- para[2]
  k1 <- para[3]
  mu2 <- para[4]
  sigma2 <- para[5]
  k2 <- para[6]
  a <- para[7]
  b <- para[8]
  res <- NULL
  try({
    res <-
      (nls(
        r ~ k1 * exp(-1 / 2 * (x - mu1) ^ 2 / sigma1 ^ 2) +
          k2 * exp(-1 / 2 * (x - mu2) ^ 2 / sigma2 ^ 2) +
          a * x + b,
        start = c(
          k1 = k1, mu1 = mu1, sigma1 = sigma1,
          k2 = k2, mu2 = mu2, sigma2 = sigma2,
          a = a, b = x[1]
        ) ,  # a = 1e-10,
        data = tab
      ))
  }, silent = T)
  if (is.null(res)) {
    cat(sprintf('nls curve fit failed @ \n'))
    return(NULL)
  }
  v <- summary(res)$parameters[,"Estimate"]
  fit <-
    function(x)
      v[1] * exp(-1 / 2 * (x - v[2]) ^ 2 / v[3] ^ 2) +
    v[4] * exp(-1 / 2 * (x - v[5]) ^ 2 / v[6] ^ 2) +
    v[7] * x + v[8]
  y <- fit(seq_along(x))
  if (!area.output)
    return(y)
  fit.baseline <- function(x)
    v[4] * x + v[5]
  baseline <- fit.baseline(x)
  area <- finiteIntegrate(x = x, y = y - baseline)
  return(area)
}