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.r...@gmail.com> wrote: > > On Wed, 1 Jul 2020 15:24:35 +0200 > Luigi Marongiu <marongiu.lu...@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 ______________________________________________ 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.