thr3ads.net - R help - [R] Problem with nls.lm function of minpack.lm package. [Mar 2011]

If this information is useful, please help other people find it:
Share via:

Yann Périard Larrivée

2011-Mar-15 15:30 UTC

[R] Problem with nls.lm function of minpack.lm package.

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
>
>
>-------------- section suivante --------------
Un texte encapsul? et encod? dans un jeu de caract?res inconnu a ?t? nettoy?...
Nom : mod.direct.txt
URL :
<https://stat.ethz.ch/pipermail/r-help/attachments/20110315/f544feb4/attachment.txt>

Ravi Varadhan

2011-Mar-15 16:28 UTC

head link

[R] Problem with nls.lm function of minpack.lm package.

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.

Maybe Matching Threads

Search for more possibly parallel threads

R help - Mar 2011 - Problem with nls.lm function of minpack.lm package.

[R] Problem with nls.lm function of minpack.lm package.

[R] Problem with nls.lm function of minpack.lm package.

Maybe Matching Threads