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.elli...@lgcgroup.com> wrote:
>
> Ivan Krylov [krylov.r...@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}}
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
