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

[Luigi Marongiu]

Got it, thanks!
Now I get:
```
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_old <- 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)
# a =  4380.4979, b = 30.3779
opt_estim <- c(4.741739,   18.905561,   42.309262,   74.655637,  115.541787,
               164.471381,  220.869196,
284.097173,  353.471198,  428.277856,  507.790487,  591.283989,
678.047947,  767.397828,
858.684086,  951.299179, 1044.682566, 1138.323858, 1231.764323,
1324.596988 , 1416.465585,
1507.062590, 1596.126574, 1683.439081, 1768.821202, 1852.130003,
1933.254918, 2012.114210,
2088.651563, 2162.832870, 2234.643232, 2304.084201, 2371.171268,
2435.931609, 2498.402055,
2558.627308, 2616.658366, 2672.551143, 2726.365284, 2778.163130,
2828.008847, 2875.967677,
2922.105311, 2966.487363, 3009.178936, 3050.244270, 3089.746448,
3127.747177, 3164.306598,
3199.483163, 3233.333529, 3265.912492, 3297.272946, 3327.465861,
3356.540283, 3384.543344,
3411.520287, 3437.514499, 3462.567553, 3486.719252)

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("bottomright",
       legend = c("Actual values", "Raw estimate", "Optimized values"),
       lty = c(1, 2, 3), lwd = c(2, 1,1))
```

And it works even better with the function Gomperz:
```
gomp = function(p, x, y) (p$a * exp(-p$b * exp(-p$c * x))) - y
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 = 3.344e+03, b = 7.715e+00, c = 1.007e-01
GP = c(3.123603,    6.093680,   11.150677,   19.256792,   31.559926,
49.332695,
       73.883838,
106.453532,  148.107305,  199.643054,  261.522677,  333.835070,
416.291985,  508.253731,
608.778470,  716.687199,  830.636345,  949.190754, 1070.891363,
1194.313689, 1318.114912,
1441.068876, 1562.089431, 1680.243298, 1794.754104, 1904.999383,
2010.502340, 2110.919990,
2206.029092, 2295.711021, 2379.936466, 2458.750617, 2532.259301,
2600.616339, 2664.012281,
2722.664577, 2776.809146, 2826.693281, 2872.569775, 2914.692156,
2953.310891, 2988.670438,
3021.007015, 3050.546986, 3077.505743, 3102.087011, 3124.482483,
3144.871725, 3163.422286,
3180.289965, 3195.619209, 3209.543577, 3222.186275, 3233.660724,
3244.071144, 3253.513139,
3262.074288, 3269.834707, 3276.867604, 3283.239800)
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))
```
Thank you

On Wed, Jul 1, 2020 at 3:42 PM Ivan Krylov <krylov.r00t using gmail.com> wrote:
>
> On Wed, 1 Jul 2020 15:24:35 +0200
> Luigi Marongiu <marongiu.luigi using gmail.com> wrote:
>
> >You are right: The vector X is actually Y -- the response I would like
> >to fit the curve upon. I understood I should fit nls.lm with a
> >function that describes the data (Holling or Gomperz), initial
> >parameters, and the actual values (Y). In return, I get the optimized
> >values for the parameters. But when I re-run the function that
> >describes the data with the optimized parameters, I get values close
> >to zero.
>
> The function to be minimized should have access to all three: the
> parameters to optimize, the independent and dependent variable values.
> Only then there is enough information to compute residuals and minimise
> their sum of squares.
>
> holly <- function(p, x, y) (p$a * x^2) / (p$b^2 + x^2) - y
> #               residuals =      y.hat(x, params)      - y.actual
> O <- nls.lm(par = list(a = 3261, b = 10), fn = holly, x = X, y = Y)
> summary(O)
> # Parameters:
> #    Estimate Std. Error t value Pr(>|t|)
> # a 4380.4979   106.8144   41.01   <2e-16 ***
> # b   30.3779     0.9995   30.39   <2e-16 ***
> # ---
> # Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> #
> # Residual standard error: 157.5 on 58 degrees of freedom
> # Number of iterations to termination: 7
> # Reason for termination: Relative error in the sum of squares is at
> # most `ftol'.
>
> In our previous examples we ended up asking R to minimise y.hat(x)
> calculated at wrong x values instead of minimising residuals.
>
> --
> Best regards,
> Ivan



-- 
Best regards,
Luigi