In finding optimal parameters in nonlinear optimization and nonlinear least squares problems, we frequently wish to fix one or more parameters while allowing the rest to be adjusted to explore or optimize an objective function.
This vignette discusses some ideas about specifying the fixed
parameters. A lot of the material is drawn from Nash J C (2014)
Nonlinear parameter optimization using R tools
Chichester UK: Wiley, in particular chapters 11 and 12. There is,
however, additional material concerning ways to manage extensible
models, as well as some update to the package nlsr. The
algorithm has been marginally altered to allow for different
sub-variants to be used, and mechanisms for specifying parameter
constraints and providing Jacobian approximations have been changed.
Here are some of the ways fixed parameters may be specified in R packages.
Function nlxb() in package nlsr has
argument masked:
This approach has a simplicity that is attractive, but introduces an
extra argument to calling sequences. (This approach was previously in
defunct package nlmrt.)
Simlarly, function nlfb() in nlsr has
argument maskidx:
From Rvmmin and Rcgmin in package
optimx the argument bdmsk:
Note that the function bmchk() in package
optimx contains a much more extensive examination of the
bounds on parameters. In particular, it considers the issues of
inadmissible bounds (lower > upper), when to convert a pair of bounds
where upper[“parameter”] - lower[“parameter”] < tol to a fixed or
masked parameter (maskadded) and whether parameters outside
of bounds should be moved to the nearest bound
(parchanged). It may be useful to use
inadmissible to refer to situations where a lower bound
is higher than an upper bound and infeasible where a
parameter value, especially in a given starting vector, is outside the
bounds.
Further in package optimx, the function
optimr() can call many different “optimizers” (actually
function minimization methods that may include bounds and possibly
masks). These may be specified by setting the lower and upper bounds
equal for the parameters to be fixed. This seems a simple method for
specifying masks, but does pose some issues. For example, what happens
when the upper bound is only very slightly greater than the lower bound?
Also should we stop or declare an error if starting values are NOT on
the fixed value?
Of these methods, my preference is now to use the last one – setting lower and upper bounds equal, and furthermore requiring the starting value of such a parameter to this fixed value, otherwise declaring an error. The approach does not add any special argument for masking, and is relatively obvious to novice users. However, such users may be tempted to put in narrow bounds rather than explicit equalities, and this could have deleterious consequences.
In the revision to package nlsr, package
nlsr, I have stopped using masked in
nlxb() and maskidx in nlfb()
(though the latter is a returned value). This is because I feel the use
of equal lower and upper bounds is a better approach. Moreover, though
it is not documented, it appears to “mostly work” for the base R
function nls() with the algorithm="port"
option and with minpack.lm::nlsLM().
bdmsk is the internal structure used in
Rcgmin, Rvmmin and nlfb to handle
bounds constraints as well as masks. There is one element of
bdmsk for each parameter, and in Rcgmin and
Rvmmin, this is used on input to specify parameter i as
fixed or masked by setting bdmsk[i] <- 0. Free
parameters have their bdmsk element 1, but
during optimization in the presence of bounds, we can set other values.
The full set is as follows
Not all these possibilities will be used by all methods that use
bdmsk.
The -1 and -3 are historical, and arose in the development of BASIC
codes for Nash and Walker-Smith (1987)
(This is now available for free download from archive.org. (https://archive.org/details/NLPE87plus). In particular,
adding 2 to the bdmsk element gives 1 for an upper bound
and -1 for a lower bound, simplifying the expression to decide if an
optimization trial step will move away from a bound.
Because masks (fixed parameters) reduce the dimensionality of the optimization problem, we can consider modifying the problem to the lower dimension space. This is Duncan Murdoch’s suggestion, using
fn0(par0) to be the initial user function of the full
dimension parameter vector par0fn1(par1) to be the reduced or internal functin of the
reduced dimension vector par1par1 <- forward(par0)par0 <- inverse(par1)The major advantage of this approach is explicit dimension reduction. The main disadvantage is the effort of transformation at every step of an optimization.
An alternative is to use the bdmsk vector to
mask the optimization search or adjustment vector,
including gradients and (approximate) Hessian or Jacobian matrices. A 0
element of bdmsk “multiplies” any adjustment. The principal
difficulty is to ensure we do not essentially divide by zero in applying
any inverse Hessian. This approach avoids forward,
inverse and fn1. However, it may hide the
reduction in dimension, and caution is necessary in using the function
and its derived gradient, Hessian and derived information.
## Loading required package: optimxsq<-function(x){
nn<-length(x)
yy<-1:nn
f<-sum((yy-x)^2)
f
}
sq.g <- function(x){
nn<-length(x)
yy<-1:nn
gg<- 2*(x - yy)
}
xx <- c(.3, 4)
uncans <- Rvmmin(xx, sq, sq.g)
proptimr(uncans)## Result uncans ( -> ) calc. min. = 0 at
## 1 2 NA NA NA NA NA NA
## After 4 fn evals, and 3 gr evals and NA hessian evals
## Termination code is 2 : Rvmminu appears to have converged
##
## -------------------------------------------------## Result cans ( -> ) calc. min. = 0.49 at
## 0.3 2 NA NA NA NA NA NA
## After 6 fn evals, and 4 gr evals and NA hessian evals
## Termination code is 2 : Rvmminb appears to have converged
##
## -------------------------------------------------## Loading required package: nlsrweed <- c(5.308, 7.24, 9.638, 12.866, 17.069, 23.192, 31.443,
38.558, 50.156, 62.948, 75.995, 91.972)
ii <- 1:12
wdf <- data.frame(weed, ii)
weedux <- nlxb(weed~b1/(1+b2*exp(-b3*ii)), start=c(b1=200, b2=50, b3=0.3))
weedux## $resid
## [1] 0.01189993 -0.03275547 0.09202995 0.20878182 0.39263404 -0.05759436
## [7] -1.10572842 0.71578576 -0.10764762 -0.34839635 0.65259251 -0.28756791
## attr(,"gradient")
## b1 b2 b3
## [1,] 0.02711658 -0.1054282 5.175642
## [2,] 0.03673674 -0.1414187 13.884948
## [3,] 0.04959588 -0.1883714 27.742382
## [4,] 0.06664474 -0.2485844 48.813666
## [5,] 0.08900539 -0.3240359 79.537273
## [6,] 0.11792062 -0.4156794 122.438293
## [7,] 0.15463505 -0.5224121 179.522463
## [8,] 0.20018622 -0.6398588 251.293721
## [9,] 0.25510631 -0.7594104 335.526319
## [10,] 0.31908250 -0.8682775 426.251654
## [11,] 0.39068787 -0.9513292 513.725387
## [12,] 0.46733360 -0.9948174 586.046595
##
## $jacobian
## b1 b2 b3
## [1,] 0.02711658 -0.1054282 5.175642
## [2,] 0.03673674 -0.1414187 13.884948
## [3,] 0.04959588 -0.1883714 27.742382
## [4,] 0.06664474 -0.2485844 48.813666
## [5,] 0.08900539 -0.3240359 79.537273
## [6,] 0.11792062 -0.4156794 122.438293
## [7,] 0.15463505 -0.5224121 179.522463
## [8,] 0.20018622 -0.6398588 251.293721
## [9,] 0.25510631 -0.7594104 335.526319
## [10,] 0.31908250 -0.8682775 426.251654
## [11,] 0.39068787 -0.9513292 513.725387
## [12,] 0.46733360 -0.9948174 586.046595
##
## $feval
## [1] 6
##
## $jeval
## [1] 6
##
## $coefficients
## b1 b2 b3
## 196.1862616 49.0916394 0.3135697
##
## $ssquares
## [1] 2.587277
##
## $lower
## [1] -Inf -Inf -Inf
##
## $upper
## [1] Inf Inf Inf
##
## $maskidx
## integer(0)
##
## $weights
## [1] 1 1 1 1 1 1 1 1 1 1 1 1
##
## $formula
## weed ~ b1/(1 + b2 * exp(-b3 * ii))
##
## $resfn
## function (prm)
## {
## if (is.null(names(prm)))
## names(prm) <- names(pvec)
## localdata <- list2env(as.list(prm), parent = data)
## eval(residexpr, envir = localdata)
## }
## <bytecode: 0x5581e0f3e468>
## <environment: 0x5581e0edcf00>
##
## $data
## parent.frame()
##
## attr(,"class")
## [1] "nlsr"## $resid
## [1] 0.01189993 -0.03275547 0.09202995 0.20878182 0.39263404 -0.05759436
## [7] -1.10572842 0.71578576 -0.10764762 -0.34839635 0.65259251 -0.28756791
## attr(,"gradient")
## b1 b2 b3
## [1,] 0.02711658 -0.1054282 5.175642
## [2,] 0.03673674 -0.1414187 13.884948
## [3,] 0.04959588 -0.1883714 27.742382
## [4,] 0.06664474 -0.2485844 48.813666
## [5,] 0.08900539 -0.3240359 79.537273
## [6,] 0.11792062 -0.4156794 122.438293
## [7,] 0.15463505 -0.5224121 179.522463
## [8,] 0.20018622 -0.6398588 251.293721
## [9,] 0.25510631 -0.7594104 335.526319
## [10,] 0.31908250 -0.8682775 426.251654
## [11,] 0.39068787 -0.9513292 513.725387
## [12,] 0.46733360 -0.9948174 586.046595
##
## $jacobian
## b1 b2 b3
## [1,] 0.02711658 -0.1054282 5.175642
## [2,] 0.03673674 -0.1414187 13.884948
## [3,] 0.04959588 -0.1883714 27.742382
## [4,] 0.06664474 -0.2485844 48.813666
## [5,] 0.08900539 -0.3240359 79.537273
## [6,] 0.11792062 -0.4156794 122.438293
## [7,] 0.15463505 -0.5224121 179.522463
## [8,] 0.20018622 -0.6398588 251.293721
## [9,] 0.25510631 -0.7594104 335.526319
## [10,] 0.31908250 -0.8682775 426.251654
## [11,] 0.39068787 -0.9513292 513.725387
## [12,] 0.46733360 -0.9948174 586.046595
##
## $feval
## [1] 6
##
## $jeval
## [1] 6
##
## $coefficients
## b1 b2 b3
## 196.1862616 49.0916394 0.3135697
##
## $ssquares
## [1] 2.587277
##
## $lower
## [1] -Inf -Inf -Inf
##
## $upper
## [1] Inf Inf Inf
##
## $maskidx
## integer(0)
##
## $weights
## [1] 1 1 1 1 1 1 1 1 1 1 1 1
##
## $formula
## weed ~ b1/(1 + b2 * exp(-b3 * ii))
##
## $resfn
## function (prm)
## {
## if (is.null(names(prm)))
## names(prm) <- names(pvec)
## localdata <- list2env(as.list(prm), parent = data)
## eval(residexpr, envir = localdata)
## }
## <bytecode: 0x5581e0f3e468>
## <environment: 0x5581dfdfb388>
##
## $data
## parent.frame()
##
## attr(,"class")
## [1] "nlsr"rfn <- function(bvec, weed=weed, ii=ii){
res <- rep(NA, length(ii))
for (i in ii){
res[i]<- bvec[1]/(1+bvec[2]*exp(-bvec[3]*i))-weed[i]
}
res
}
weeduf <- nlfb(start=c(200, 50, 0.3),resfn=rfn,weed=weed, ii=ii,
control=list(japprox="jacentral"))
weeduf## $resid
## [1] 0.01189992 -0.03275548 0.09202995 0.20878181 0.39263404 -0.05759436
## [7] -1.10572842 0.71578576 -0.10764762 -0.34839635 0.65259251 -0.28756791
##
## $jacobian
## [,1] [,2] [,3]
## [1,] 0.02711658 -0.1054282 5.175642
## [2,] 0.03673674 -0.1414187 13.884948
## [3,] 0.04959588 -0.1883714 27.742382
## [4,] 0.06664474 -0.2485844 48.813666
## [5,] 0.08900539 -0.3240359 79.537273
## [6,] 0.11792062 -0.4156794 122.438293
## [7,] 0.15463505 -0.5224121 179.522463
## [8,] 0.20018622 -0.6398588 251.293721
## [9,] 0.25510631 -0.7594104 335.526319
## [10,] 0.31908250 -0.8682775 426.251654
## [11,] 0.39068787 -0.9513292 513.725387
## [12,] 0.46733360 -0.9948174 586.046595
## attr(,"gradient")
## [,1] [,2] [,3]
## [1,] 0.02711658 -0.1054282 5.175642
## [2,] 0.03673674 -0.1414187 13.884948
## [3,] 0.04959588 -0.1883714 27.742382
## [4,] 0.06664474 -0.2485844 48.813666
## [5,] 0.08900539 -0.3240359 79.537273
## [6,] 0.11792062 -0.4156794 122.438293
## [7,] 0.15463505 -0.5224121 179.522463
## [8,] 0.20018622 -0.6398588 251.293721
## [9,] 0.25510631 -0.7594104 335.526319
## [10,] 0.31908250 -0.8682775 426.251654
## [11,] 0.39068787 -0.9513292 513.725387
## [12,] 0.46733360 -0.9948174 586.046595
##
## $feval
## [1] 6
##
## $jeval
## [1] 6
##
## $coefficients
## p1 p2 p3
## 196.1862616 49.0916394 0.3135697
##
## $ssquares
## [1] 2.587277
##
## $lower
## [1] -Inf -Inf -Inf
##
## $upper
## [1] Inf Inf Inf
##
## $maskidx
## integer(0)
##
## $weights0
## NULL
##
## $weights
## [1] 1 1 1 1 1 1 1 1 1 1 1 1
##
## $formula
## NULL
##
## attr(,"class")
## [1] "nlsr"weedcf <- nlfb(start=c(200, 50, 0.3),resfn=rfn,weed=weed, ii=ii, lower=c(200, 0, 0),
upper=c(200, 100,100), control=list(japprox="jacentral"))
weedcf## $resid
## [1] 0.05960669 0.01784926 0.14230513 0.25376273 0.42559950 -0.04444931
## [7] -1.11938060 0.67287349 -0.17220893 -0.41045953 0.63860337 -0.18656556
##
## $jacobian
## [,1] [,2] [,3]
## [1,] 0 -0.1055032 5.223551
## [2,] 0 -0.1412715 13.988935
## [3,] 0 -0.1878788 27.906100
## [4,] 0 -0.2476049 49.036487
## [5,] 0 -0.3224405 79.821472
## [6,] 0 -0.4134148 122.811031
## [7,] 0 -0.5196038 180.082070
## [8,] 0 -0.6369429 252.284531
## [9,] 0 -0.7572465 337.427048
## [10,] 0 -0.8681500 429.828206
## [11,] 0 -0.9547420 519.970633
## [12,] 0 -1.0030657 595.951255
## attr(,"gradient")
## [,1] [,2] [,3]
## [1,] 0 -0.1055032 5.223551
## [2,] 0 -0.1412715 13.988935
## [3,] 0 -0.1878788 27.906100
## [4,] 0 -0.2476049 49.036487
## [5,] 0 -0.3224405 79.821472
## [6,] 0 -0.4134148 122.811031
## [7,] 0 -0.5196038 180.082070
## [8,] 0 -0.6369429 252.284531
## [9,] 0 -0.7572465 337.427048
## [10,] 0 -0.8681500 429.828206
## [11,] 0 -0.9547420 519.970633
## [12,] 0 -1.0030657 595.951255
##
## $feval
## [1] 4
##
## $jeval
## [1] 4
##
## $coefficients
## p1 p2 p3
## 200.0000000 49.5108198 0.3114607
##
## $ssquares
## [1] 2.618154
##
## $lower
## [1] 200 0 0
##
## $upper
## [1] 200 100 100
##
## $maskidx
## [1] 1
##
## $weights0
## NULL
##
## $weights
## [1] 1 1 1 1 1 1 1 1 1 1 1 1
##
## $formula
## NULL
##
## attr(,"class")
## [1] "nlsr"Package nlraa has a selfStart model SSbell
(Archontoulis and Miguez (2013)) of which
the formula is
y ≈ ymax * exp(a * (x − xc)2 + b * (x − xc)3)
This is essentially the Gaussian bell curve with an additional cubic element in the exponential function. If we fix b = 0, then we have the usual Gaussian, and we can use the standard deviation sigma of the variable x with xc equal to its mean and our parameters will be given approximately by
ymax = max(y) a = −0.5/sigma2 xc = mean(y)
We illustrate this in the following example.
## Loading required package: ggplot2set.seed(1234)
x <- 1:20
y <- bell(x, 8, -0.0314, 0.000317, 13) + rnorm(length(x), 0, 0.5)
dat <- data.frame(x = x, y = y)
fit <- nls(y ~ SSbell(x, ymax, a, b, xc), data = dat)
fit## Nonlinear regression model
## model: y ~ SSbell(x, ymax, a, b, xc)
## data: dat
## ymax a b xc
## 7.7601430 -0.0311315 0.0005088 13.0872979
## residual sum-of-squares: 4.517
##
## Number of iterations to convergence: 6
## Achieved convergence tolerance: 5.664e-06gfrm<- y ~ ymax * exp(a *(x - xc)^2)
gfrmx<- y ~ ymax * exp(a *(x - xc)^2 + b*(x - xc)^3)
stgauss<-c(ymax=max(y), a=-0.5/(sd(x)^2), xc=mean(x))
cat("stgauss:"); print(stgauss)## stgauss:## ymax a xc
## 7.78739029 -0.01428571 10.50000000st2<-c(ymax=8, a= 0.03, xc= 13)
st3<-c(ymax=8, a= 0.03, b=0, xc= 13)
fit2 <- nls(gfrm, start=st2, data=dat, trace=TRUE)## 493716.4 (4.02e+01): par = (8 0.03 13)
## 6154.601 (4.38e+00): par = (0.9881319 0.0300116 12.82135)
## 1234.014 (1.98e+00): par = (0.9909919 0.02862245 11.63345)
## 311.0878 (1.50e+00): par = (1.748681 0.01208336 11.58072)
## 101.1911 (3.03e+00): par = (4.421808 -0.01330662 10.76144)
## 66.60083 (3.15e+00): par = (6.385106 -0.02233438 15.9676)
## 20.36513 (1.83e+00): par = (6.317481 -0.02153419 12.63178)
## 5.558036 (4.18e-01): par = (7.606842 -0.03007893 13.50289)
## 4.716550 (2.95e-02): par = (7.765594 -0.03131867 13.25154)
## 4.712500 (8.56e-04): par = (7.789061 -0.03157114 13.25752)
## 4.712497 (2.95e-05): par = (7.789119 -0.03157216 13.257)
## 4.712497 (1.13e-06): par = (7.789137 -0.03157246 13.257)##
## Formula: y ~ ymax * exp(a * (x - xc)^2)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## ymax 7.789137 0.246937 31.54 < 2e-16 ***
## a -0.031572 0.002471 -12.78 3.84e-10 ***
## xc 13.256996 0.146867 90.27 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5265 on 17 degrees of freedom
##
## Number of iterations to convergence: 11
## Achieved convergence tolerance: 1.133e-06# Use the selfStart model
xfrm <- y ~ SSbell(x, ymax, a, b, xc)
library(nlsr)
fx2mx <- nlxb(xfrm, start=st3, data=dat, lower=c(0, -1e5, 0, 0),
upper=c(1e5, 1e5, 0, 1e5), trace=TRUE, control=list(japprox="jacentral"))## The following parameters are masked:[1] "b"
## Using approximation jacentral
## is.character(ctrl$japprox)= TRUE
## No backtrack
## The following parameters are masked:[1] "b"
## 00lamda: 1e-04 SS= 493716.4 ( NA ) at ymax = 8 a = 0.03 b = 0 xc = 13 f/j 1 / 0
## <<lamda: 4e-05 SS= 6154.66 ( 0.9995642 ) at ymax = 0.9881376 a = 0.0300116 b = 0 xc = 12.82135 f/j 2 / 1
## <<lamda: 1.6e-05 SS= 1234.022 ( 0.9644597 ) at ymax = 0.9909919 a = 0.02862244 b = 0 xc = 11.63346 f/j 3 / 2
## <<lamda: 6.4e-06 SS= 311.0898 ( 0.8142677 ) at ymax = 1.748675 a = 0.01208346 b = 0 xc = 11.58072 f/j 4 / 3
## <<lamda: 2.56e-06 SS= 101.1908 ( 0.8091139 ) at ymax = 4.42179 a = -0.01330648 b = 0 xc = 10.76147 f/j 5 / 4
## <<lamda: 1.024e-06 SS= 66.60145 ( 0.7980163 ) at ymax = 6.385122 a = -0.02233452 b = 0 xc = 15.96761 f/j 6 / 5
## <<lamda: 4.096e-07 SS= 20.36515 ( 0.760554 ) at ymax = 6.317463 a = -0.02153393 b = 0 xc = 12.63181 f/j 7 / 6
## <<lamda: 1.6384e-07 SS= 5.558011 ( 0.5599207 ) at ymax = 7.606844 a = -0.03007892 b = 0 xc = 13.50289 f/j 8 / 7
## <<lamda: 6.5536e-08 SS= 4.71655 ( 0.3355534 ) at ymax = 7.765595 a = -0.03131868 b = 0 xc = 13.25154 f/j 9 / 8
## <<lamda: 2.62144e-08 SS= 4.7125 ( 0.02349175 ) at ymax = 7.789061 a = -0.03157114 b = 0 xc = 13.25752 f/j 10 / 9
## <<lamda: 1.048576e-08 SS= 4.712497 ( 0.0007767177 ) at ymax = 7.789119 a = -0.03157216 b = 0 xc = 13.257 f/j 11 / 10
## <<lamda: 4.194304e-09 SS= 4.712497 ( 2.632395e-05 ) at ymax = 7.789137 a = -0.03157246 b = 0 xc = 13.257 f/j 12 / 11
## <<lamda: 1.677722e-09 SS= 4.712497 ( 9.104232e-07 ) at ymax = 7.789137 a = -0.03157247 b = 0 xc = 13.257 f/j 13 / 12## $resid
## [1] 0.62108737 -0.11360768 -0.51326828 1.19439823 -0.22124270 -0.31655543
## [7] 0.13804413 0.02112411 -0.06570564 0.03814849 -0.16732284 0.15905068
## [13] 0.16103846 -0.13283967 -0.47668090 0.11489604 0.32048064 0.48764713
## [19] 0.40173622 -1.26892833
## attr(,"gradient")
## ymax a b xc
## [1,] 0.008710237 10.1926565 -124.9313480 -0.05250998
## [2,] 0.018299866 18.0626827 -203.3315430 -0.10132072
## [3,] 0.036094609 29.5782199 -303.3836774 -0.18209179
## [4,] 0.066836434 44.6110812 -412.9645916 -0.30430647
## [5,] 0.116187850 61.7012815 -509.4672230 -0.47185725
## [6,] 0.189620226 77.7835474 -564.4748776 -0.67681408
## [7,] 0.290526070 88.5943811 -554.3346709 -0.89408503
## [8,] 0.417890080 89.9552810 -472.8945351 -1.08050690
## [9,] 0.564307105 79.6546814 -339.0896446 -1.18153501
## [10,] 0.715394305 59.1111242 -192.5246837 -1.14601560
## [11,] 0.851435870 33.7833559 -76.2488927 -0.94517135
## [12,] 0.951338145 11.7082479 -14.7172186 -0.58816150
## [13,] 0.997916911 0.5133763 -0.1319356 -0.12613868
## [14,] 0.982721274 4.2257347 3.1397386 0.35912817
## [15,] 0.908537457 21.4995338 37.4737775 0.77887743
## [16,] 0.788554692 46.2140579 126.7653547 1.06386408
## [17,] 0.642535765 70.1176376 262.4506116 1.18289303
## [18,] 0.491517853 86.1262734 408.4972757 1.14662302
## [19,] 0.352986252 90.6828999 520.7922744 0.99706800
## [20,] 0.237986756 84.2847476 568.3324065 0.78928540
##
## $jacobian
## [,1] [,2] [,3] [,4]
## [1,] 0.008710237 10.1926565 0 -0.05250998
## [2,] 0.018299866 18.0626827 0 -0.10132072
## [3,] 0.036094609 29.5782199 0 -0.18209179
## [4,] 0.066836434 44.6110812 0 -0.30430647
## [5,] 0.116187850 61.7012815 0 -0.47185725
## [6,] 0.189620226 77.7835473 0 -0.67681408
## [7,] 0.290526070 88.5943811 0 -0.89408503
## [8,] 0.417890080 89.9552810 0 -1.08050690
## [9,] 0.564307105 79.6546814 0 -1.18153501
## [10,] 0.715394305 59.1111241 0 -1.14601560
## [11,] 0.851435871 33.7833559 0 -0.94517135
## [12,] 0.951338146 11.7082478 0 -0.58816150
## [13,] 0.997916912 0.5133763 0 -0.12613868
## [14,] 0.982721275 4.2257347 0 0.35912817
## [15,] 0.908537458 21.4995337 0 0.77887743
## [16,] 0.788554692 46.2140579 0 1.06386408
## [17,] 0.642535765 70.1176376 0 1.18289303
## [18,] 0.491517854 86.1262733 0 1.14662302
## [19,] 0.352986253 90.6828998 0 0.99706800
## [20,] 0.237986756 84.2847476 0 0.78928540
## attr(,"gradient")
## [,1] [,2] [,3] [,4]
## [1,] 0.008710237 10.1926565 0 -0.05250998
## [2,] 0.018299866 18.0626827 0 -0.10132072
## [3,] 0.036094609 29.5782199 0 -0.18209179
## [4,] 0.066836434 44.6110812 0 -0.30430647
## [5,] 0.116187850 61.7012815 0 -0.47185725
## [6,] 0.189620226 77.7835473 0 -0.67681408
## [7,] 0.290526070 88.5943811 0 -0.89408503
## [8,] 0.417890080 89.9552810 0 -1.08050690
## [9,] 0.564307105 79.6546814 0 -1.18153501
## [10,] 0.715394305 59.1111241 0 -1.14601560
## [11,] 0.851435871 33.7833559 0 -0.94517135
## [12,] 0.951338146 11.7082478 0 -0.58816150
## [13,] 0.997916912 0.5133763 0 -0.12613868
## [14,] 0.982721275 4.2257347 0 0.35912817
## [15,] 0.908537458 21.4995337 0 0.77887743
## [16,] 0.788554692 46.2140579 0 1.06386408
## [17,] 0.642535765 70.1176376 0 1.18289303
## [18,] 0.491517854 86.1262733 0 1.14662302
## [19,] 0.352986253 90.6828998 0 0.99706800
## [20,] 0.237986756 84.2847476 0 0.78928540
##
## $feval
## [1] 13
##
## $jeval
## [1] 13
##
## $coefficients
## ymax a b xc
## 7.78913698 -0.03157247 0.00000000 13.25699581
##
## $ssquares
## [1] 4.712497
##
## $lower
## [1] 0e+00 -1e+05 0e+00 0e+00
##
## $upper
## [1] 1e+05 1e+05 0e+00 1e+05
##
## $maskidx
## [1] 3
##
## $weights
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## $formula
## y ~ SSbell(x, ymax, a, b, xc)
##
## $resfn
## function (prm)
## {
## if (is.null(names(prm)))
## names(prm) <- names(pvec)
## localdata <- list2env(as.list(prm), parent = data)
## eval(residexpr, envir = localdata)
## }
## <bytecode: 0x5581e0f3e468>
## <environment: 0x5581ebe35ae8>
##
## $data
## dat
##
## attr(,"class")
## [1] "nlsr"# Or the formula for the same model with Jacobian approximation
fx2mgc <- nlxb(gfrmx, start=st3, data=dat, lower=c(0, -1e5, 0, 0),
upper=c(1e5, 1e5, 0, 1e5), trace=TRUE, control=list(japprox="jacentral"))## The following parameters are masked:[1] "b"
## Using approximation jacentral
## is.character(ctrl$japprox)= TRUE
## No backtrack
## The following parameters are masked:[1] "b"
## 00lamda: 1e-04 SS= 493716.4 ( NA ) at ymax = 8 a = 0.03 b = 0 xc = 13 f/j 1 / 0
## <<lamda: 4e-05 SS= 6154.66 ( 0.9995642 ) at ymax = 0.9881376 a = 0.0300116 b = 0 xc = 12.82135 f/j 2 / 1
## <<lamda: 1.6e-05 SS= 1234.022 ( 0.9644597 ) at ymax = 0.9909919 a = 0.02862244 b = 0 xc = 11.63346 f/j 3 / 2
## <<lamda: 6.4e-06 SS= 311.0898 ( 0.8142677 ) at ymax = 1.748675 a = 0.01208346 b = 0 xc = 11.58072 f/j 4 / 3
## <<lamda: 2.56e-06 SS= 101.1908 ( 0.8091139 ) at ymax = 4.42179 a = -0.01330648 b = 0 xc = 10.76147 f/j 5 / 4
## <<lamda: 1.024e-06 SS= 66.60145 ( 0.7980163 ) at ymax = 6.385122 a = -0.02233452 b = 0 xc = 15.96761 f/j 6 / 5
## <<lamda: 4.096e-07 SS= 20.36515 ( 0.760554 ) at ymax = 6.317463 a = -0.02153393 b = 0 xc = 12.63181 f/j 7 / 6
## <<lamda: 1.6384e-07 SS= 5.558011 ( 0.5599207 ) at ymax = 7.606844 a = -0.03007892 b = 0 xc = 13.50289 f/j 8 / 7
## <<lamda: 6.5536e-08 SS= 4.71655 ( 0.3355534 ) at ymax = 7.765595 a = -0.03131868 b = 0 xc = 13.25154 f/j 9 / 8
## <<lamda: 2.62144e-08 SS= 4.7125 ( 0.02349175 ) at ymax = 7.789061 a = -0.03157114 b = 0 xc = 13.25752 f/j 10 / 9
## <<lamda: 1.048576e-08 SS= 4.712497 ( 0.0007767177 ) at ymax = 7.789119 a = -0.03157216 b = 0 xc = 13.257 f/j 11 / 10
## <<lamda: 4.194304e-09 SS= 4.712497 ( 2.632395e-05 ) at ymax = 7.789137 a = -0.03157246 b = 0 xc = 13.257 f/j 12 / 11
## <<lamda: 1.677722e-09 SS= 4.712497 ( 9.104232e-07 ) at ymax = 7.789137 a = -0.03157247 b = 0 xc = 13.257 f/j 13 / 12## $resid
## [1] 0.62108737 -0.11360768 -0.51326828 1.19439823 -0.22124270 -0.31655543
## [7] 0.13804413 0.02112411 -0.06570564 0.03814849 -0.16732284 0.15905068
## [13] 0.16103846 -0.13283967 -0.47668090 0.11489604 0.32048064 0.48764713
## [19] 0.40173622 -1.26892833
##
## $jacobian
## [,1] [,2] [,3] [,4]
## [1,] 0.008710237 10.1926565 0 -0.05250998
## [2,] 0.018299866 18.0626827 0 -0.10132072
## [3,] 0.036094609 29.5782199 0 -0.18209179
## [4,] 0.066836434 44.6110812 0 -0.30430647
## [5,] 0.116187850 61.7012815 0 -0.47185725
## [6,] 0.189620226 77.7835473 0 -0.67681408
## [7,] 0.290526070 88.5943811 0 -0.89408503
## [8,] 0.417890080 89.9552810 0 -1.08050690
## [9,] 0.564307105 79.6546814 0 -1.18153501
## [10,] 0.715394305 59.1111241 0 -1.14601560
## [11,] 0.851435871 33.7833559 0 -0.94517135
## [12,] 0.951338146 11.7082478 0 -0.58816150
## [13,] 0.997916912 0.5133763 0 -0.12613868
## [14,] 0.982721275 4.2257347 0 0.35912817
## [15,] 0.908537458 21.4995337 0 0.77887743
## [16,] 0.788554692 46.2140579 0 1.06386408
## [17,] 0.642535765 70.1176376 0 1.18289303
## [18,] 0.491517854 86.1262733 0 1.14662302
## [19,] 0.352986253 90.6828998 0 0.99706800
## [20,] 0.237986756 84.2847476 0 0.78928540
## attr(,"gradient")
## [,1] [,2] [,3] [,4]
## [1,] 0.008710237 10.1926565 0 -0.05250998
## [2,] 0.018299866 18.0626827 0 -0.10132072
## [3,] 0.036094609 29.5782199 0 -0.18209179
## [4,] 0.066836434 44.6110812 0 -0.30430647
## [5,] 0.116187850 61.7012815 0 -0.47185725
## [6,] 0.189620226 77.7835473 0 -0.67681408
## [7,] 0.290526070 88.5943811 0 -0.89408503
## [8,] 0.417890080 89.9552810 0 -1.08050690
## [9,] 0.564307105 79.6546814 0 -1.18153501
## [10,] 0.715394305 59.1111241 0 -1.14601560
## [11,] 0.851435871 33.7833559 0 -0.94517135
## [12,] 0.951338146 11.7082478 0 -0.58816150
## [13,] 0.997916912 0.5133763 0 -0.12613868
## [14,] 0.982721275 4.2257347 0 0.35912817
## [15,] 0.908537458 21.4995337 0 0.77887743
## [16,] 0.788554692 46.2140579 0 1.06386408
## [17,] 0.642535765 70.1176376 0 1.18289303
## [18,] 0.491517854 86.1262733 0 1.14662302
## [19,] 0.352986253 90.6828998 0 0.99706800
## [20,] 0.237986756 84.2847476 0 0.78928540
##
## $feval
## [1] 13
##
## $jeval
## [1] 13
##
## $coefficients
## ymax a b xc
## 7.78913698 -0.03157247 0.00000000 13.25699581
##
## $ssquares
## [1] 4.712497
##
## $lower
## [1] 0e+00 -1e+05 0e+00 0e+00
##
## $upper
## [1] 1e+05 1e+05 0e+00 1e+05
##
## $maskidx
## [1] 3
##
## $weights
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## $formula
## y ~ ymax * exp(a * (x - xc)^2 + b * (x - xc)^3)
##
## $resfn
## function (prm)
## {
## if (is.null(names(prm)))
## names(prm) <- names(pvec)
## localdata <- list2env(as.list(prm), parent = data)
## eval(residexpr, envir = localdata)
## }
## <bytecode: 0x5581e0f3e468>
## <environment: 0x5581ec2c3c50>
##
## $data
## dat
##
## attr(,"class")
## [1] "nlsr"# Or the formula and analytic derivatives
fx2mg <- nlxb(gfrmx, start=st3, data=dat, lower=c(0, -1e5, 0, 0),
upper=c(1e5, 1e5, 0, 1e5), trace=TRUE)## The following parameters are masked:[1] "b"
## No backtrack
## The following parameters are masked:[1] "b"
## 00lamda: 1e-04 SS= 493716.4 ( NA ) at ymax = 8 a = 0.03 b = 0 xc = 13 f/j 1 / 0
## <<lamda: 4e-05 SS= 6154.66 ( 0.9995642 ) at ymax = 0.9881376 a = 0.0300116 b = 0 xc = 12.82135 f/j 2 / 1
## <<lamda: 1.6e-05 SS= 1234.022 ( 0.9644597 ) at ymax = 0.9909919 a = 0.02862244 b = 0 xc = 11.63346 f/j 3 / 2
## <<lamda: 6.4e-06 SS= 311.0898 ( 0.8142677 ) at ymax = 1.748675 a = 0.01208346 b = 0 xc = 11.58072 f/j 4 / 3
## <<lamda: 2.56e-06 SS= 101.1908 ( 0.8091139 ) at ymax = 4.42179 a = -0.01330648 b = 0 xc = 10.76147 f/j 5 / 4
## <<lamda: 1.024e-06 SS= 66.60145 ( 0.7980163 ) at ymax = 6.385122 a = -0.02233452 b = 0 xc = 15.96761 f/j 6 / 5
## <<lamda: 4.096e-07 SS= 20.36515 ( 0.760554 ) at ymax = 6.317463 a = -0.02153393 b = 0 xc = 12.63181 f/j 7 / 6
## <<lamda: 1.6384e-07 SS= 5.558011 ( 0.5599207 ) at ymax = 7.606844 a = -0.03007892 b = 0 xc = 13.50289 f/j 8 / 7
## <<lamda: 6.5536e-08 SS= 4.71655 ( 0.3355534 ) at ymax = 7.765595 a = -0.03131868 b = 0 xc = 13.25154 f/j 9 / 8
## <<lamda: 2.62144e-08 SS= 4.7125 ( 0.02349175 ) at ymax = 7.789061 a = -0.03157114 b = 0 xc = 13.25752 f/j 10 / 9
## <<lamda: 1.048576e-08 SS= 4.712497 ( 0.0007767176 ) at ymax = 7.789119 a = -0.03157216 b = 0 xc = 13.257 f/j 11 / 10
## <<lamda: 4.194304e-09 SS= 4.712497 ( 2.632404e-05 ) at ymax = 7.789137 a = -0.03157246 b = 0 xc = 13.257 f/j 12 / 11
## <<lamda: 1.677722e-09 SS= 4.712497 ( 9.104065e-07 ) at ymax = 7.789137 a = -0.03157247 b = 0 xc = 13.257 f/j 13 / 12## $resid
## [1] 0.62108737 -0.11360768 -0.51326828 1.19439823 -0.22124270 -0.31655543
## [7] 0.13804413 0.02112411 -0.06570564 0.03814849 -0.16732284 0.15905068
## [13] 0.16103846 -0.13283967 -0.47668090 0.11489604 0.32048064 0.48764713
## [19] 0.40173622 -1.26892834
## attr(,"gradient")
## ymax a b xc
## [1,] 0.008710237 10.1926565 -124.9313479 -0.05250998
## [2,] 0.018299866 18.0626827 -203.3315430 -0.10132072
## [3,] 0.036094609 29.5782199 -303.3836773 -0.18209179
## [4,] 0.066836434 44.6110812 -412.9645915 -0.30430647
## [5,] 0.116187850 61.7012815 -509.4672230 -0.47185725
## [6,] 0.189620226 77.7835474 -564.4748775 -0.67681408
## [7,] 0.290526070 88.5943811 -554.3346709 -0.89408503
## [8,] 0.417890080 89.9552810 -472.8945351 -1.08050690
## [9,] 0.564307105 79.6546814 -339.0896446 -1.18153501
## [10,] 0.715394305 59.1111242 -192.5246837 -1.14601560
## [11,] 0.851435870 33.7833559 -76.2488927 -0.94517135
## [12,] 0.951338145 11.7082479 -14.7172186 -0.58816150
## [13,] 0.997916911 0.5133763 -0.1319356 -0.12613868
## [14,] 0.982721274 4.2257347 3.1397386 0.35912817
## [15,] 0.908537457 21.4995338 37.4737775 0.77887743
## [16,] 0.788554692 46.2140579 126.7653547 1.06386408
## [17,] 0.642535765 70.1176376 262.4506116 1.18289303
## [18,] 0.491517853 86.1262734 408.4972757 1.14662302
## [19,] 0.352986252 90.6828999 520.7922744 0.99706800
## [20,] 0.237986756 84.2847476 568.3324065 0.78928540
##
## $jacobian
## ymax a b xc
## [1,] 0.008710237 10.1926565 0 -0.05250998
## [2,] 0.018299866 18.0626827 0 -0.10132072
## [3,] 0.036094609 29.5782199 0 -0.18209179
## [4,] 0.066836434 44.6110812 0 -0.30430647
## [5,] 0.116187850 61.7012815 0 -0.47185725
## [6,] 0.189620226 77.7835474 0 -0.67681408
## [7,] 0.290526070 88.5943811 0 -0.89408503
## [8,] 0.417890080 89.9552810 0 -1.08050690
## [9,] 0.564307105 79.6546814 0 -1.18153501
## [10,] 0.715394305 59.1111242 0 -1.14601560
## [11,] 0.851435870 33.7833559 0 -0.94517135
## [12,] 0.951338145 11.7082479 0 -0.58816150
## [13,] 0.997916911 0.5133763 0 -0.12613868
## [14,] 0.982721274 4.2257347 0 0.35912817
## [15,] 0.908537457 21.4995338 0 0.77887743
## [16,] 0.788554692 46.2140579 0 1.06386408
## [17,] 0.642535765 70.1176376 0 1.18289303
## [18,] 0.491517853 86.1262734 0 1.14662302
## [19,] 0.352986252 90.6828999 0 0.99706800
## [20,] 0.237986756 84.2847476 0 0.78928540
##
## $feval
## [1] 13
##
## $jeval
## [1] 13
##
## $coefficients
## ymax a b xc
## 7.78913698 -0.03157247 0.00000000 13.25699581
##
## $ssquares
## [1] 4.712497
##
## $lower
## [1] 0e+00 -1e+05 0e+00 0e+00
##
## $upper
## [1] 1e+05 1e+05 0e+00 1e+05
##
## $maskidx
## [1] 3
##
## $weights
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## $formula
## y ~ ymax * exp(a * (x - xc)^2 + b * (x - xc)^3)
##
## $resfn
## function (prm)
## {
## if (is.null(names(prm)))
## names(prm) <- names(pvec)
## localdata <- list2env(as.list(prm), parent = data)
## eval(residexpr, envir = localdata)
## }
## <bytecode: 0x5581e0f3e468>
## <environment: 0x5581ec5a89b0>
##
## $data
## dat
##
## attr(,"class")
## [1] "nlsr"## fx2mx -- ss= 4.712497 : ymax = 7.789137 a = -0.03157247 b = 0 xc = 13.257; 13 res/ 13 jac## fx2mgc -- ss= 4.712497 : ymax = 7.789137 a = -0.03157247 b = 0 xc = 13.257; 13 res/ 13 jac## fx2mg -- ss= 4.712497 : ymax = 7.789137 a = -0.03157247 b = 0 xc = 13.257; 13 res/ 13 jac