Search code examples
pythonrgaussian

Is there a function for a self start model to estimate the parameter of a Gaussian model in R?


I need to estimate the parameters of a Gaussian model. I like to use self start to select the parameters. I already did something like this using a logistic model with a different dataset. I would like to know if there is a similar function of SSlogis for gaussian

x <- runif(20)
y <- rexp(20)
a <- data.frame(x,y)

log_model <- nls(y~SSlogis(x, phi1, phi2, phi3), data = a)

What I want is to do a similar approach but using a gaussian model. I didn't find a function similar to SSlogis for a gaussian distribution. The R documentation said that there is an SSgauss function but I didn't find it https://www.rdocumentation.org/packages/xcms/versions/1.48.0/topics/SSgauss I saw something similar to what I want in a python3 code.

def gaussian_f(x,a,b,c):
   y = a * np.exp(-0.5 * ((x-b)/c)**2)
   return y
##optimize from scipy
gaussian_m, cov = optimize.curve_fit(gaussian_f, x=np.arrange(len(a["y"])), y=dtf["y"].values, maxfev=10000, p0=[1,np.mean(a["y"]) ,1]

Solution

  • R is open source and Bioconductor releases package xcms under

    License: GPL (>= 2)

    so the source code for function SSgauss can be downloaded and used as long as the users follow the licensing terms.

    Here is the source code for the function found in file xcms/R/models.R.

    SSgauss <- selfStart(~ h*exp(-(x-mu)^2/(2*sigma^2)), function(mCall, data, LHS) {
      
      xy <- sortedXyData(mCall[["x"]], LHS, data)
      
      len <- dim(xy)[1]
      xyarea <- sum((xy[2:len,2]+xy[1:(len-1),2])*(xy[2:len,1]-xy[1:(len-1),1]))/2
      maxpos <- which.max(xy[,2])
      
      mu <- xy[maxpos,1]
      h <- xy[maxpos,2]
      sigma <- xyarea/(h*sqrt(2*pi))
      
      value <- c(mu, sigma, h)
      names(value) <- mCall[c("mu", "sigma", "h")]
      value
      
    }, c("mu", "sigma", "h"))
    

    Now fit a gaussian model to the data set in the question.

    x <- runif(20)
    y <- rexp(20)
    a <- data.frame(x, y)
    
    gauss_model <- nls(y ~ SSgauss(x, mu, sigma, h), data = a)
    summary(gauss_model)
    #
    #Formula: y ~ SSgauss(x, mu, sigma, h)
    #
    #Parameters:
    #      Estimate Std. Error t value Pr(>|t|)    
    #mu      0.5844     0.0989   5.909 1.72e-05 ***
    #sigma   0.3540     0.1436   2.465  0.02463 *  
    #h       1.2453     0.3364   3.702  0.00177 ** 
    #---
    #Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    #
    #Residual standard error: 0.7832 on 17 degrees of freedom
    #
    #Number of iterations to convergence: 9 
    #Achieved convergence tolerance: 2.897e-06