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