[R] argument "x" is missing in minpack.lm

[Luigi Marongiu]

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 using gmail.com> wrote:
>
> (Adding R-help back to Cc:)
>
> On Tue, 30 Jun 2020 14:44:29 +0200
> Luigi Marongiu <marongiu.luigi using 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