One things you should do, as Kate suggested, is to check whether the Jacobian functiions are correctly code. You can do this with "numDeriv" package:
require(numDeriv) ?jacobian Compare the jacobian from numDeriv with your jacobian for a few reasonable parameter vectors. Ravi. ____________________________________________________________________ Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvarad...@jhmi.edu ----- Original Message ----- From: Yann Périard Larrivée <yann.periard-larrive...@ulaval.ca> Date: Tuesday, March 15, 2011 11:57 am Subject: [R] Problem with nls.lm function of minpack.lm package. To: "r-help@r-project.org" <r-help@r-project.org> > Dear R useRs, > > I have a problem with nls.lm function of minpackl.lm package. > > I need to fit the Van Genuchten Model to a set of data of Theta and > hydraulic conductivity with nls.lm function of minpack.lm package. > > For the first fit, the parameter estimates keep changing even after > 1000 iterations (Th) > > and > > I have a following error message for fit of hydraulic conductivity (k); > > Reason for termination: The cosine of the angle between `fvec' and > any column of the > Jacobian is at most `gtol' in absolute value. > > I dont know what I can do for solve my problem > > Yesterday I asked to Katharine Mullen > > Here are her aswer > > Yann Périard > > Étudiant à la maîtrise en sols et environnement > Département des sols et de génie agroalimentaire > > Centre de Recherche en Horticulture > Pavillon Envirotron, Université Laval > 2480 boulevard Hochelaga, local 2241 > Québec (Qc), Canada > G1V 0A6 > Tél.: 418-656-2131 poste 8229 > Courriel : yann.periard-larrive...@ulaval.ca > ________________________________________ > De : Katharine Mullen [kmul...@nist.gov] > Date d'envoi : 14 mars 2011 18:05 > À : Yann Périard Larrivée > Objet : Re: your mail > > Dear Yann Périard, > > Thanks for including a reproducible example. The error message you see > generally means that nls.lm cannot improve the model fit further, and > so > is stopping. This happens when the objective function being minimized > (the sum squared residuals) does not change as the algorithm tries to > vary > the parameters (the way the parameter vector changes is decided by either > a finite difference method or the Jacobian function you supply). > > There are some strange things going on with your models; in the case > of > this call: > > out <- nls.lm(par = guess, fn = fcn, jac = fcn.jac, > fcall = f.Th, jcall = j.Th, > h = h, Th = Th, control = nls.lm.control(maxfev = integer(), > maxiter = 10000, nprint=100)) > > The parameter estimates keep changing even after 1000 iterations. > This is > not good. > > In the second fit, I change the starting values slightly, to > > ##starting values (alp = 0.04, n = 1.6, L = 0.5) > guess.k <- c(alp = 0.08, n = 1.61, L = 0.51) > > ## to use an analytical expression for the gradient found in fcn.jac > ## uncomment jac = fcn.jac > out.k <- nls.lm(par = guess.k, fn = fcn.k, jac = fcn.jac.k, > fcall = f.k, jcall = j.k, > h = h, k = k, control = nls.lm.control(gtol = 1, > epsfcn = > .01, > factor = 100, maxiter = 10000, nprint > = 50)) > > and do not get any movement of the parameters toward the correct values. > > I suggest that you try to figure out what is wrong with your Jacobian > functions and/or model functions. Also consider posting the message > you > sent me to the R-help mailing list; many people that read this list are > interested in nonlinear regression problems, and will be eager to offer > their insights. > > Regarding combination of two residual functions: using finite difference > derivatives, it is straightforward. You define a new residual function, > resid2, that evaluates both old residual functions, and returns their > scaled sums. So if you had > > foo1(par1, ...) > > and foo2(par2, ...) > > you define > > resid2 <- function(c(par1,par2), ...) { > sse1 <- foo1(par1, ...) > sse2 <- foo2(par2, ...) > return(sse1 * w1 + sse2 * w2) > } > > w1 and w2 are weights that determine how much each of the old residual > functions contributes to the combined fit. > > Hope this helps, and best regards, > Kate > > -- > Katharine Mullen > Structure Determination Methods Group, Ceramics Division > National Institute of Standards and Technology (NIST) > 100 Bureau Drive, M/S 8520 > Gaithersburg, MD, 20899, USA > phone: 001 301 975 6890 > email: katharine.mul...@nist.gov > > On Mon, 14 Mar 2011, Yann Périard Larrivée wrote: > > > Hi Mrs Mullen > > > > I have a problem for fitting the following Van Genuchten Model to a > set of data of > > hydraulic conductivity with nls.lm function of minpack.lm package. > > I have a following error message; > > > > Reason for termination: The cosine of the angle between `fvec' and > any column of the > > Jacobian is at most `gtol' in absolute value. > > I dont know what I can do for solve my problem > > > > Here his the code > > > > #fonction pour simulation des données de conductivité hydraulique > > f.k <- function(h, alp, n, L, ksat){ > > ksat <- 0.000108 > > k.mod > <-expression(ksat+((((1+(alp*h)^n)^(n-(1/n))-((alp*h)^(n-1)))^2)/((1+(alp*h)^n)^((n-(1/n))*(L+2))))) > > > > eval (k) > > } > > #fonction d'aide pour le gradient analytique > > j.k <- function(h, alp, n, L, ksat){ > > ksat <- 0.000108 > > k.mod > <-expression(ksat+((((1+(alp*h)^n)^(n-(1/n))-((alp*h)^(n-1)))^2)/((1+(alp*h)^n)^((n-(1/n))*(L+2))))) > > > > c(eval(D(k.mod, "alp")), eval(D(k.mod, "n")), eval(D(k.mod, "L"))) > > } > > h <- hydro$h > > # Paramètre surestimé pour démarer la simulation > > p.k <- c(alp = 0.04, n = 1.6, L = 0.5) > > ## get data with noise > > k <- hydro$k > > ## plot the data to fit > > par(mfrow=c(2,1), mar = c(3,5,2,1)) > > plot(h, k, bg = "black", cex = 0.5, main="data", > ylab=expression(K(h)), log="y") > > ## define a residual function > > fcn.k <- function(p.k, h, k, fcall, jcall) > > (k - do.call("fcall", c(list(h = h), as.list(p.k)))) > > ##define analytical expression for the gradient > > fcn.jac.k <- function(p.k, h, k, fcall, jcall) > > -do.call("jcall", c(list(h = h), as.list(p.k))) > > ##starting values (alp = 0.04, n = 1.6, L = 0.5) > > guess.k <- c(alp = 0.04, n = 1.6, L = 0.5) > > ## to use an analytical expression for the gradient found in fcn.jac > > ## uncomment jac = fcn.jac > > out.k <- nls.lm(par = guess.k, fn = fcn.k, jac = fcn.jac.k, > > fcall = f.k, jcall = j.k, > > h = h, k = k, control = nls.lm.control(gtol = 1, > epsfcn = 0, factor = 100, > > maxiter = 100, nprint = 50)) > > ## get the fitted values > > N1.k <- do.call("f.k", c(list(h = h), out.k$par)) > > ## add a blue line representing the fitting values to the plot of data > > lines(h, N1.k, col="blue", lwd=2) > > summary(out.k) > > > > I have another question? > > > > It is possible to use two models with the same parameters (alp, n) > for fitting Th~h and k~h > > simultaneously with following models > > I guess, but I dont know how I can do that > > > > f <- function(h, Th, k, Ths, Thr, alp, n, L, ksat){ > > Th_mod <- expression(Thr+(Ths-Thr)/(1+(alp*h)^(n-1/n))) > > eval (Th_mod) > > k_mod <- > > > expression(ksat+(((1+(alp*h)^n)-((alp*h)^(n-1)))^2)/((1+(alp*h)^n)^((n-1/n)*(L+2)))) > > eval (k_mod) > > } > > #fonction d'aide pour le gradient analytique > > j <- function(h, Th, K, Ths, Thr, alp, n, L, ksat){ > > Th_mod <- expression(Thr+(Ths-Thr)/(1+(alp*h)^(n-1/n))) > > k_mod <- > > > expression(ksat+(((1+(alp*h)^n)-((alp*h)^(n-1)))^2)/((1+(alp*h)^n)^((n-1/n)*(L+2)))) > > c(eval(D(Th_mod, "Thr")), eval(D(Th_mod, k_mod, "alp")), > eval(D(Th_mod, k_mod, "n")), > > eval(D(k_mod, "L"))) > > } > > > > Thanks for your help > > > > Regards > > > > Yann Périard > > > > > > Étudiant à la maîtrise en sols et environnement > > > > Département des sols et de génie agroalimentaire > > > > > > > > Centre de Recherche en Horticulture > > > > Pavillon Envirotron, Université Laval > > > > 2480 boulevard Hochelaga, local 2241 > > > > Québec (Qc), Canada > > > > G1V 0A6 > > > > Tél.: 418-656-2131 poste 8229 > > > > Courriel : yann.periard-larrive...@ulaval.ca > > > > > > ______________________________________________ > R-help@r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.