But I get the same error even if I run:
```
parms = list(a=3261, b=10, x=CH$Cum_Dead[1:60])
O = nls.lm(parms, holling,
lower=NULL, upper=NULL, jac = NULL)> Error in fn(par, ...) : argument "x" is missing, with no default
```
the function to be optimized is:
```
function(a, b, x) {
y = (a * x^2) / (b^2 + x^2)
return(y)
}
```
so the parameters a, b, x are the same I called within nls.lm
On Tue, Jun 30, 2020 at 1:32 PM Sarah Goslee <sarah.goslee at gmail.com>
wrote:>
> You create objects A, B, and K, but then use A, B, and X in the call to
nls.lm().
>
> I have no idea if nls.lm is the right tool for your job, but I do know that
you need to use the same names.
>
> Sarah
>
> On Tue, Jun 30, 2020 at 6:01 AM Luigi Marongiu <marongiu.luigi at
gmail.com> wrote:
>>
>> Hello,
>> I am trying to optimize a function with the function nls.lm from the
>> package minpack.lm. But I can't get it right:
>> ```
>> A = 3261
>> B = 10
>> K = 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)
>> O = nls.lm(list(a=A, b=B, x=X), holling,
>> lower=NULL, upper=NULL, jac = NULL)
>> > Error in fn(par, ...) : argument "x" is missing, with no
default
>> ```
>> What am I missing?
>> Is nls.lm the right function for functions that are not linear?
>> Thank you
>> --
>> Best regards,
>> Luigi
>>
>> ______________________________________________
>> R-help at 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.
>
> --
> Sarah Goslee (she/her)
> http://www.sarahgoslee.com
--
Best regards,
Luigi
On Tue, 30 Jun 2020 13:53:10 +0200 Luigi Marongiu <marongiu.luigi at gmail.com> wrote:> function(a, b, x) { > y = (a * x^2) / (b^2 + x^2) > return(y) > }Take a look at the examples in ?nls.lm. The first argument of the function must be the parameter list/vector; the rest of the arguments are passed from ... in nls.lm call. This is _unlike_ do.call, which can be used to "unpack" a list of arguments into function arguments. To make it more concrete, change your function to: function(par) (par$a * par$x^2) / (par$b^2 + par$x^2) Then it should work. -- Best regards, Ivan
(Adding R-help back to Cc:) On Tue, 30 Jun 2020 14:44:29 +0200 Luigi Marongiu <marongiu.luigi at gmail.com> wrote:> Ok, I tried with: > ``` > holly <- function(p) { > y = (p$a * p$x^2) / (p$b^2 + p$x^2) > return(y) > } > X = 1:60 > A = 3261 > B = 10> 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)You are correct in returning a vector of residuals to minimise a sum of squares of, but X seems to be an independent variable, not a parameter to optimize, so it shouldn't be passed as such. 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 or pass X as an argument that nls.lm will pass to your function: holly <- function(p, X) (p$a * X^2) / (p$b^2 + X^2) # nls.lm will pass the X argument to function holly O <- nls.lm(par = list(a = 3261, b = 10), fn = holly, X = X) summary(O) # Parameters: # Estimate Std. Error t value Pr(>|t|) # a 3.090e-16 4.102e-17 7.533e+00 3.72e-10 *** # b 1.000e+01 1.525e-08 6.558e+08 < 2e-16 *** # --- # Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 # # Residual standard error: 3.107e-16 on 58 degrees of freedom # Number of iterations to termination: 2 # Reason for termination: Relative error between `par' and the solution # is at most `ptol'. -- Best regards, Ivan
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 XThat 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 use...{{dropped:8}}
Thank you,
I got this:
```
holly = function(p, x) {
y = (p$a * x^2) / (p$b^2 + x^2)
return(y)
}
A = 3261
B = 10
K = CH$Cum_Dead[1:60]
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)
O <- nls.lm(par = list(a = A, b = B), fn = holly, x =
X)> summary(O)
Parameters:
Estimate Std. Error t value Pr(>|t|)
a 3.090e-16 4.102e-17 7.533e+00 3.72e-10 ***
b 1.000e+01 1.525e-08 6.558e+08 < 2e-16 ***
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
Residual standard error: 3.107e-16 on 58 degrees of freedom
Number of iterations to termination: 2
Reason for termination: Relative error between `par' and the solution
is at most `ptol'.
```
Is this correct?
if yes, then the case is closed, thank you.
However, the optimization is worse the the eyeball estimate:
```
Addendum:
the optimization actually got a worse outcome than the original
eyeball estimation:
```
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 = 10
raw_estim <- c(32.28713, 125.42308, 269.25688, 449.79310, 652.20000,
863.20588, 1072.40940, 1272.58537,
1459.34254, 1630.50000, 1785.43439, 1924.52459, 2048.73234,
2159.31081, 2257.61538, 2344.98876,
2422.69666, 2491.89623, 2553.62473, 2608.80000, 2658.22736,
2702.60959, 2742.55803, 2778.60355,
2811.20690, 2840.76804, 2867.63450, 2892.10860, 2914.45377,
2934.90000, 2953.64844, 2970.87544,
2986.73591, 3001.36624, 3014.88679, 3027.40401, 3039.01225,
3049.79534, 3059.82788, 3069.17647,
3077.90062, 3086.05365, 3093.68343, 3100.83301, 3107.54118,
3113.84296, 3119.77003, 3125.35108,
3130.61216, 3135.57692, 3140.26694, 3144.70185, 3148.89962,
3152.87666, 3156.64800, 3160.22744,
3163.62765, 3166.86028, 3169.93605, 3172.86486)
# a = 3.090e-16, b = 1.000e+01
opt_estim <- c(3.059406e-18, 1.188462e-17, 2.551376e-17, 4.262069e-17,
6.180000e-17,
8.179412e-17, 1.016174e-16,
1.205854e-16, 1.382818e-16, 1.545000e-16, 1.691810e-16, 1.823607e-16,
1.941301e-16, 2.046081e-16,
2.139231e-16, 2.222022e-16, 2.295656e-16, 2.361226e-16, 2.419718e-16,
2.472000e-16, 2.518835e-16,
2.560890e-16, 2.598744e-16, 2.632899e-16, 2.663793e-16, 2.691804e-16,
2.717262e-16, 2.740452e-16,
2.761626e-16, 2.781000e-16, 2.798765e-16, 2.815089e-16, 2.830118e-16,
2.843981e-16, 2.856792e-16,
2.868653e-16, 2.879653e-16, 2.889870e-16, 2.899377e-16, 2.908235e-16,
2.916502e-16, 2.924227e-16,
2.931457e-16, 2.938232e-16, 2.944588e-16, 2.950560e-16, 2.956176e-16,
2.961464e-16, 2.966449e-16,
2.971154e-16, 2.975598e-16, 2.979800e-16, 2.983778e-16, 2.987546e-16,
2.991120e-16, 2.994512e-16,
2.997734e-16, 3.000797e-16, 3.003711e-16, 3.006486e-16)
plot(1:60, actual, lty = 1 , type = "l", lwd = 2,
xlab = "Index", ylab = "Values")
points(1:60, raw_estim, lty = 2, type = "l")
points(1:60, opt_estim, lty = 3, type = "l")
legend("right",
legend = c("Actual values", "Raw estimate",
"Optimized values"),
lty = c(1, 2, 3), lwd = c(2, 1,1))
```
Is that normal?
On Tue, Jun 30, 2020 at 3:41 PM Ivan Krylov <krylov.r00t at gmail.com>
wrote:>
> (Adding R-help back to Cc:)
>
> On Tue, 30 Jun 2020 14:44:29 +0200
> Luigi Marongiu <marongiu.luigi at gmail.com> wrote:
>
> > Ok, I tried with:
> > ```
> > holly <- function(p) {
> > y = (p$a * p$x^2) / (p$b^2 + p$x^2)
> > return(y)
> > }
> > X = 1:60
> > A = 3261
> > B = 10
>
> > 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)
>
> You are correct in returning a vector of residuals to minimise a sum of
> squares of, but X seems to be an independent variable, not a parameter
> to optimize, so it shouldn't be passed as such. 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
>
> or pass X as an argument that nls.lm will pass to your function:
>
> holly <- function(p, X) (p$a * X^2) / (p$b^2 + X^2)
> # nls.lm will pass the X argument to function holly
> O <- nls.lm(par = list(a = 3261, b = 10), fn = holly, X = X)
> summary(O)
> # Parameters:
> # Estimate Std. Error t value Pr(>|t|)
> # a 3.090e-16 4.102e-17 7.533e+00 3.72e-10 ***
> # b 1.000e+01 1.525e-08 6.558e+08 < 2e-16 ***
> # ---
> # Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
> #
> # Residual standard error: 3.107e-16 on 58 degrees of freedom
> # Number of iterations to termination: 2
> # Reason for termination: Relative error between `par' and the solution
> # is at most `ptol'.
>
> --
> Best regards,
> Ivan
--
Best regards,
Luigi