[R] Problem with nls.lm function of minpack.lm package.
Yann Périard Larrivée
yann.periard-larrivee.1 at ulaval.ca
Tue Mar 15 16:30:56 CET 2011
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>
More information about the R-help
mailing list