Addendum.
I have found the function Gompertz even better than the Holling III
because it gives more pronounced S profile. However the optimization
is bad even in this case:
```
gomp = function(p, x) {
y = p$a * exp(-p$b * exp(-p$c * x))
return(y)
}
A = 3261
B = 10
C = 1
X = c(8, 24, 39, 63, 89, 115, 153, 196, 242, 287, 344, 408, 473,
546, 619, 705, 794, 891, 999, 1096, 1242, 1363, 1506, 1648, 1753,
1851, 1987, 2101, 2219, 2328, 2425, 2575, 2646, 2698, 2727, 2771, 2818,
2853, 2895, 2926, 2964, 2995, 3025, 3053, 3080, 3102, 3119, 3141, 3152,
3159, 3172, 3182, 3196, 3209, 3220, 3231, 3239, 3246, 3252, 3261)
W = nls.lm(par = list(a = A, b = B, c = C), fn = gomp, x = X)
W
Nonlinear regression via the Levenberg-Marquardt algorithm
parameter estimates: 2.81739569520133e-17, 20.0000002056654, 0.100000000717244 residual sum-of-squares: 4.556e-32 reason terminated: Relative error between `par' and the solution is at most `ptol'. summary(W)
Parameters:
Estimate Std. Error t value Pr(>|t|)
a 2.817e-17 3.764e-18 7.485e+00 4.94e-10 ***
b 2.000e+01 1.872e-08 1.068e+09 < 2e-16 ***
c 1.000e-01 2.924e-11 3.420e+09 < 2e-16 ***
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
Residual standard error: 2.827e-17 on 57 degrees of freedom
Number of iterations to termination: 2
Reason for termination: Relative error between `par' and the solution
is at most `ptol'.
```
but eyeballing gives:
```
actual <- c(8, 24, 39, 63, 89, 115, 153, 196, 242, 287,
344, 408, 473,
546, 619, 705, 794, 891, 999, 1096, 1242, 1363,
1506, 1648, 1753,
1851, 1987, 2101, 2219, 2328, 2425, 2575, 2646,
2698, 2727, 2771, 2818,
2853, 2895, 2926, 2964, 2995, 3025, 3053, 3080,
3102, 3119, 3141, 3152,
3159, 3172, 3182, 3196, 3209, 3220, 3231, 3239,
3246, 3252, 3261)
# a = 3261, b = 20, c = 0.1
GOMP = c(4.508508e-05, 2.523166e-04, 1.198631e-03, 4.909360e-03, 1.758298e-02,
5.577422e-02, 1.585133e-01,
4.078768e-01, 9.592495e-01, 2.079652e+00, 4.188604e+00,
7.892437e+00, 1.400134e+01, 2.351997e+01,
3.760684e+01, 5.750425e+01, 8.444654e+01, 1.195592e+02,
1.637624e+02, 2.176925e+02, 2.816487e+02,
3.555714e+02, 4.390496e+02, 5.313549e+02, 6.314945e+02,
7.382759e+02, 8.503763e+02, 9.664098e+02,
1.084989e+03, 1.204775e+03, 1.324520e+03, 1.443095e+03,
1.559508e+03, 1.672916e+03, 1.782624e+03,
1.888078e+03, 1.988863e+03, 2.084686e+03, 2.175363e+03,
2.260805e+03, 2.341005e+03, 2.416022e+03,
2.485969e+03, 2.551004e+03, 2.611315e+03, 2.667115e+03,
2.718631e+03, 2.766102e+03, 2.809769e+03,
2.849874e+03, 2.886657e+03, 2.920347e+03, 2.951171e+03,
2.979341e+03, 3.005063e+03, 3.028528e+03,
3.049918e+03, 3.069402e+03, 3.087140e+03, 3.103278e+03)
# a = 2.817e-17, b = 2.000e+01, c = 1.000e-01
GP = c(3.894654e-25, 2.179626e-24, 1.035432e-23, 4.240928e-23, 1.518898e-22,
4.818030e-22, 1.369310e-21,
3.523425e-21, 8.286434e-21, 1.796498e-20, 3.618306e-20,
6.817846e-20, 1.209500e-19, 2.031762e-19,
3.248650e-19, 4.967478e-19, 7.294876e-19, 1.032806e-18,
1.414654e-18, 1.880527e-18, 2.433010e-18,
3.071587e-18, 3.792710e-18, 4.590085e-18, 5.455136e-18,
6.377563e-18, 7.345937e-18, 8.348287e-18,
9.372625e-18, 1.040739e-17, 1.144180e-17, 1.246611e-17,
1.347174e-17, 1.445141e-17, 1.539911e-17,
1.631008e-17, 1.718071e-17, 1.800847e-17, 1.879178e-17,
1.952986e-17, 2.022267e-17, 2.087070e-17,
2.147493e-17, 2.203673e-17, 2.255773e-17, 2.303975e-17,
2.348477e-17, 2.389484e-17, 2.427206e-17,
2.461851e-17, 2.493625e-17, 2.522729e-17, 2.549356e-17,
2.573691e-17, 2.595910e-17, 2.616180e-17,
2.634657e-17, 2.651489e-17, 2.666812e-17, 2.680752e-17)
plot(1:60, actual, lty = 1 , type = "l", lwd = 2,
xlab = "Index", ylab = "Values")
points(1:60, GOMP, lty = 2, type = "l")
points(1:60, GP, lty = 3, type = "l")
legend("right",
legend = c("Actual values", "Raw estimate", "Optimized values"),
lty = c(1, 2, 3), lwd = c(2, 1,1))
```
On Tue, Jun 30, 2020 at 9:59 PM Stephen Ellison <S.Ellison at lgcgroup.com> wrote:
Ivan Krylov [krylov.r00t at gmail.com] said:
Instead you can either close over X: X <- c(...) holly <- function(p) (p$a * X^2) / (p:b^2 + X^2) # function holly now "contains" the vector X
That would not be an accurate statement as written.
The function only contains an unevaluated call referencing X; not the vector X.
If X is not defined inside the function or its arguments, scoping rules take over and R goes looking in the function's environment, using the first thing called "X" that it finds.
So
X<-1:5
h <- function(p=2) p*X
h()
# works, but
X <- pi
h()
# Not the same answer. If the function 'contained' the vector, the result would be 2*(1:5) as above.
# This is why it's not wise to rely on scoping rules in functions, unless you _completely_ control the environment.
# and
rm(X)
h()
# returns an error because X is no longer defined, in the function or out of it, even though X was defined at the time h() was defined.
BUT
X <- pi/2
fh <- function(p=2) {
X <- 7
h(p)
}
fh()
# returns pi and not 14 because h() was bound to the global environment on creation and still is when R finds it on evaluating h() in the fh() function body.
Moral: if you want to be safe, pass it as an argument.
*******************************************************************
This email and any attachments are confidential. Any u...{{dropped:14}}