myprojects-hg-reborn

Fork

rm(list=ls())
library(minpack.lm)
set.seed(123) # I set the seed of the rng in order to have reproducible results
f <- function(T, tau, N0, a, f0) {
    expr <- expression(N0*exp(-T/tau)*(1 + a*cos(f0*T)))
    eval(expr) ### expression I want to use for my NLS
}
j <- function(T, tau, N0, a, f0) {
    expr <- expression(N0*exp(-T/tau)*(1 + a*cos(f0*T)))
    c(eval(D(expr, "tau")),   #the call eval(D(...)) means that R has to evaluate symbo
                              #lically the derivatives of the function (this is done here to
                              # get the jacobian)
      eval(D(expr, "N0" )),
      eval(D(expr, "a" )),      # expression evaluated with the trial parameters
      eval(D(expr, "f0" )))
}


T <- seq(0, 8, len=501) #I generate a time sequence
p <- c("tau" = 2.2, "N0" = 1000, "a" = 0.25, "f0" = 8) #vector with the "true"
                                                       #parameters of the system
N <- do.call("f", c(list(T = T), as.list(p))) # This means that I am calling the
                                              # function f with argument and parameters given 
                                              # by T & p  

 N <- rnorm(length(N), mean=N, sd=sqrt(N)) # this line drops some Gaussian noise on the 
                                           # damped oscillations. I am not redefining N ---
                                           # subtleties of the do.call command

plot(T, N, bg = "black", pch = 21, cex = 0.5)  # I use time T as the x coord
                                              # whereas N is along the y-axis
                                              ### this is basically a plot of the data
                                              # I want to fit.         

fcn     <- function(p, T, N, N.Err, fcall, jcall)
    (N - do.call("fcall", c(list(T = T), as.list(p))))/N.Err  # Now I define the function I want
                                                        # to minimize as the difference between 
                                                        # the simulated data and the function
                                  #I plane to use to model them. 


fcn.jac <- function(p, T, N, N.Err, fcall, jcall) {
    N.Err <- rep(N.Err, length(p))                            #I need this function if I also
    -do.call("jcall", c(list(T = T), as.list(p)))/N.Err   # want to use the jacobian in the 
}                                                         # minimization (unused in the example)

guess <- c("tau" = 2.2, "N0" = 1500, "a" = 0.25, "f0" = 10) # vector with the initial guesses
                                                           # of the fitting parameters


out <- nls.lm(par = guess, fn = fcn, #jac = fcn.jac,
              fcall = f,  jcall = j,
              T = T, N = N, N.Err = sqrt(N),
              control = list(nprint = 3, diag = numeric()))

#The previous call is important. the first argument is a list of the initial values for the
# fitting parameters. the 2nd one is the function I need to minimize (difference between fcall
# and the experimental data. the third one tells the expression of fcall. The third one specifies jcall which is needed if I want to use the jacobian in the optimization
# the other arguments are clear; N.Err is chosen as sqrt(N)).
#In this case, even setting N.Err=1 has a small effect on the final results of the optimization
# since I am not using the jacobian, I can comment the part "jcall=j" without affecting the
#optimization



 N1 <- do.call("f", c(list(T = T), out$par)) # N1 == N - sqrt(N)*out$fvec
lines(T, N1, col="blue", lwd=2)
str(out)
print(sqrt(diag(solve(out$hessian)))) # calculating of SSE
# the following line is a comment
#rm(f,j,fcn,fcn.jac,T,p,guess,N,N1,out)