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