library(minpack.lm)

Non-Linear Curve Fit Example from nls.lm package

Use Levenberg-Marquardt algorithm for non-linear curve fitting.

Created simulated data with noise

Function to simulate data

f <- function(T, tau, N0, a, f0) 
{
  expr <- expression(N0 * exp(-T/tau) * (1 + a*cos(f0*T)))
  eval(expr)
}

Values over which to simulate data: 500 intervals between 0 and 8

T <- seq(0, 8, length=501)

Parameter values underlying simulated data

p <- c(tau = 2.2, N0 = 1000, a = 0.25, f0 = 8)

Define vector with simulated data to fit

Ndet <- do.call("f", c(list(T = T), as.list(p)))

Add noise. Set random number seed so example is reproducible.

set.seed(17)
N <- Ndet +  rnorm(length(Ndet), mean=Ndet, sd=.01*max(Ndet))
head(N)
## [1] 2487.312 2476.828 2444.928 2400.202 2375.693 2313.235

Non-linear curve fit

Helper function for an analytical gradient.

Function D computes derivatives with respective to named parameter.

\[ N_0 exp(-T/tau) (1 + a\cos(f0 T)) \]

See info about Hessian/Jacobian matrix.

j <- function(T, tau, N0, a, f0) 
{
  expr <- expression(N0*exp(-T/tau)*(1 + a*cos(f0*T)))
  c(eval(D(expr, "tau")), 
    eval(D(expr, "N0" )),
    eval(D(expr, "a"  )), 
    eval(D(expr, "f0" )))
}

Define a residual function

fcn <- function(p, T, N, fcall, jcall)
    (N - do.call("fcall", c(list(T = T), as.list(p))))

Define analytical expression for the gradient

fcn.jac <- function(p, T, N, fcall, jcall) 
{
  -do.call("jcall", c(list(T = T), as.list(p)))
}

Initial guess starting values

guess <- c(tau = 2.2, N0 = 1500, a = 0.25, f0 = 10)

Non-linear curve fit

out <- nls.lm(par = guess, fn = fcn, jac = fcn.jac,
              fcall = f, jcall = j,
              T = T, N = N, control = nls.lm.control(nprint=1))
## It.    0, RSS = 2.91785e+07, Par. =        2.2       1500       0.25         10
## It.    1, RSS = 9.4718e+06, Par. =     2.1019    2041.65 -0.0272045    9.91903
## It.    2, RSS = 8.88376e+06, Par. =    2.16873    2022.07  0.0487361    10.6175
## It.    3, RSS = 8.79912e+06, Par. =     2.1579    2028.62   0.028564    10.4468
## It.    4, RSS = 8.7256e+06, Par. =    2.15996    2027.57  0.0320865    10.1236
## It.    5, RSS = 8.5294e+06, Par. =    2.15985    2027.34  0.0399979    9.70825
## It.    6, RSS = 7.80537e+06, Par. =    2.15871    2027.19  0.0566179    9.13837
## It.    7, RSS = 4.03412e+06, Par. =    2.15784    2025.81  0.0996119    8.30636
## It.    8, RSS = 2.10099e+06, Par. =    2.18278    2008.66   0.225579    7.64731
## It.    9, RSS =     257598, Par. =    2.16802    2015.77   0.221037    7.95009
## It.   10, RSS =    86804.1, Par. =    2.19638    2001.55   0.251485    7.99715
## It.   11, RSS =    86020.1, Par. =    2.19702    2001.44   0.251983    7.99114
## It.   12, RSS =    86020.1, Par. =    2.19703    2001.43   0.251994    7.99114
## It.   13, RSS =    86020.1, Par. =    2.19703    2001.43   0.251994    7.99114

get the fitted values

N1 <- do.call("f", c(list(T = T), out$par))

Plots

Plot the original data with fitted solution, along with residuals.

plot(T, N, cex = 0.5, main="Simulated Data to Fit")
grid()
lines(T, N1, col="blue", lwd=2)   # Fitted solution

Add a plot of the log residual sum of squares as it is made to decrease each iteration; note that the RSS at the starting parameter values is also stored

plot(1:(out$niter+1), log(out$rsstrace), type="b",
     main="log residual sum of squares vs. iteration number",
     xlab="iteration", 
     ylab="log residual sum of squares", pch=21) 
grid()

Information regarding standard errors

summary(out) 

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
tau 2.197e+00  3.447e-03   637.3   <2e-16 ***
N0  2.001e+03  2.209e+00   906.0   <2e-16 ***
a   2.520e-01  1.137e-03   221.6   <2e-16 ***
f0  7.991e+00  2.902e-03  2754.0   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.16 on 497 degrees of freedom
Number of iterations to termination: 13 
Reason for termination: Conditions for `info = 1' and `info = 2' both hold. 
packageVersion("minpack.lm")
## [1] '1.1.8'

efg

2015-03-08 1419