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