myprojects-hg-reborn

Fork

rm(list=ls())

library(minpack.lm)


data <- read.table("number_molecules_vs_time.dat", header=FALSE)
data <- as.matrix(data)

# now exponential fitting of the maxima evolution

N<-data[ ,2]
x <- data[ ,1]

A1 <- data[1 ,2]


#I use a simple exponential function (decay) for my fittings

f <- function(x, A1, lambda) {
    expr <- expression(
A1*exp(-x/lambda)
) 
    eval(expr) ### expression I want to use for my NLS
}

#I now introduce the Jacobian of the previous function

## j <- function(x, A1, lambda) {
##     expr <- expression(
## A1*exp(-x/lambda)
##          )
##     c(eval(D(expr, "A1")),
## eval(D(expr, "lambda" ))
## )
## }


#In this case I have only a fit parameter, namely the decay rate

j <- function(x, A1, lambda) {
    expr <- expression(
A1*exp(-x/lambda)
         )
    c(eval(D(expr, "lambda" )))
}


#the function above contains the Jacobian of the function I want to use to carry out my fittings

fcn     <- function(p, x, N, N.Err, fcall, jcall)
    (N - do.call("fcall", c(list(x = x, A1=A1), 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 (eventually). 

fcn.jac <- function(p, x, N, N.Err, fcall, jcall) {
    N.Err <- rep(N.Err, length(p))                   
    -do.call("jcall", c(list(x = x, A1=A1), as.list(p)))/N.Err   
}



 myN.Err<-1



guess <- c("lambda"=3)



out <- nls.lm(par = guess, fn = fcn,   jac = fcn.jac,
              fcall = f,   jcall = j,
              x = x , N = N, N.Err = myN.Err,
              control = list(nprint = 3,maxfew=200,factor=100,ftol=1e-10, diag = numeric()))
            #   control = list(maxfew=200,factor=100,ftol=0.00001, diag = numeric()))



# now I plot some results


newx<- x
M1 <- do.call("f", c(list(x = newx, A1=A1), out$par))


pdf("absorption_fitting.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)
plot(x, N, col="blue", pch = 5,lty=1, cex.axis=1.4,cex.lab = 1.6,xlab=expression("Time"),
ylab=expression("Number of tracked molecules")
,ylim=range(c(0,max(N))))
lines(newx, M1, col="red", lwd=2.)
dev.off()





print("So far so good")