R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("minpack.lm-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('minpack.lm')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "nls.lm"
> 
> ### * nls.lm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: nls.lm
> ### Title: Solving NLS problem by Levenberg-Marquardt algorithm
> ### Aliases: nls.lm
> ### Keywords: nonlinear optimize
> 
> ### ** Examples
> 
> f <- function(T, tau, N0, a, f0) {
+     expr <- expression(N0*exp(-T/tau)*(1 + a*cos(f0*T)))
+     eval(expr)
+ }
> 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" )))
+ }
> 
> T <- seq(0, 8, len=501)
> p <- c("tau" = 2.2, "N0" = 1000, "a" = 0.25, "f0" = 8)
> 
> N <- do.call("f", c(list(T = T), as.list(p)))
> N <- rnorm(length(N), mean=N, sd=sqrt(N))
> 
> plot(T, N, bg = "black", pch = 21, cex = 0.5)
> 
> fcn     <- function(p, T, N, N.Err, fcall, jcall)
+     (N - do.call("fcall", c(list(T = T), as.list(p))))/N.Err
> fcn.jac <- function(p, T, N, N.Err, fcall, jcall) {
+     N.Err <- rep(N.Err, length(p))
+     -do.call("jcall", c(list(T = T), as.list(p)))/N.Err
+ }
> 
> guess <- c("tau" = 2.2, "N0" = 1500, "a" = 0.25, "f0" = 10)
> 
> 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()))
Iter =   0        2.2       1500       0.25         10
Iter =   3     2.1767    977.627  0.0144358    9.66853
Iter =   6     2.1736    979.558  0.0360898    8.91798
Iter =   9    2.17474    1002.07   0.173726    7.59207
Iter =  12    2.18181    1005.73   0.247692    7.99103
Iter =  14    2.18197    1005.69   0.250978    7.99242
> 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)
List of 6
 $ par    : Named num [1:4]    2.182 1005.694    0.251    7.992
  ..- attr(*, "names")= chr [1:4] "tau" "N0" "a" "f0"
 $ hessian: num [1:4, 1:4] 41062.9    55.1   107.2   199.8    55.1 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:4] "tau" "N0" "a" "f0"
  .. ..$ : chr [1:4] "tau" "N0" "a" "f0"
 $ fvec   : num [1:501] -0.8633 -0.0422 -1.0707  1.3466  0.1155 ...
 $ info   : int 1
 $ message: chr "Relative error in the sum of squares is at most `ftol'."
 $ diag   : Named num [1:4] 310.493   0.378 396.058 261.829
  ..- attr(*, "names")= chr [1:4] "tau" "N0" "a" "f0"
> 
> print(sqrt(diag(solve(out$hessian))))  # calculating of SSE
        tau          N0           a          f0 
0.007406514 4.116781060 0.003764164 0.005811234 
> #rm(f,j,fcn,fcn.jac,T,p,guess,N,N1,out)
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  0.11 0.02 0.14 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')