________________________________________
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.