R/fit.norm.R

In yanxianUCSB/yxhelper: Yx's Personal R packages

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

    # Verify para
    if (class(para) != "numeric" | length(para) <= 4) {
        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
    mu <- para[1]
    sigma <- para[2]
    k <- para[3]
    a <- para[4]
    b <- para[5]
    res <- NULL
    try({
        res <-
            (nls(
                r ~ k * exp(-1 / 2 * (x - mu) ^ 2 / sigma ^ 2) + a * x + b,  # +a*x
                start = c(
                    mu = mu, sigma = sigma, k = k, 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[3] * exp(-1 / 2 * (x - v[1]) ^ 2 / v[2] ^ 2) +  v[4] * x + v[5]  # v[4]*x +
    y <- fit(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)
}