[R] Problem with nls.lm function of minpack.lm package.
Ravi Varadhan
rvaradhan at jhmi.edu
Tue Mar 15 17:28:01 CET 2011
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: rvaradhan at jhmi.edu
----- Original Message -----
From: Yann Périard Larrivée <yann.periard-larrivee.1 at ulaval.ca>
Date: Tuesday, March 15, 2011 11:57 am
Subject: [R] Problem with nls.lm function of minpack.lm package.
To: "r-help at r-project.org" <r-help at 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-larrivee.1 at ulaval.ca
> ________________________________________
> De : Katharine Mullen [kmullen at 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.Mullen at 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-larrivee.1 at ulaval.ca
> >
> > >
> ______________________________________________
> R-help at r-project.org mailing list
>
> PLEASE do read the posting guide
> and provide commented, minimal, self-contained, reproducible code.
More information about the R-help
mailing list