I have been working on a curve fitting script which fits 3 exponentially modified Gaussians (EMGs) to a convolved curve. My base function is similar to a Gaussian distribution, but includes a third parameter (the first two being mu
and sigma
) which determines the weight of the exponential component of the function.
So overall, each EMG peak takes 3 parameters, plus an amplitude coefficient (in order to match experimental data with values > 1.0)
To deconvolve 3 EMG peaks the total # of parameters to be minimized is 3x4 = 12
In some cases the fit works nicely, but in many cases it fails to converge, and returns an error like this
Convergence failure: false convergence (8)
after only 50 or so iterations (well below the limit).
By using the trace option, I can see that the result is very close to matching the data. And by plotting the curve from my initial estimate values, it can also be seen that the starting parameters are within reasonable proximity to the data:
Data = black (noise added), initial = orange, final iteration before error = red
Here is a sample of my code, where I call nls()
:
t <- 0.05
fit <- nls(y ~ emgmix(a,b,c,d,a1,b1,c1,d1,a2,b2,c2,d2),
start = list(
a=pk1coef[2],
b=pk1coef[2],
c=t,
d=y[pk1c[2]]*40,
a1=pk2coef[1],
b1=pk2coef[2],
c1=t,
d1=y[pk2c[2]]*40,
a2=pk3coef[1],
b2=pk3coef[2],
c2=t,
d2=y[pk3c[2]]*40),
lower=rep(0.001,12),
control = list(maxiter = 1000),
trace = TRUE,
algorithm = "port",
)
so why am I getting an error when the algorithm seems to be succeeding?
0: 562831.45: 341.700 10.6000 0.0500000 27623.1 419.300 10.8000 0.0500000 2132.38 497.000 14.1000 0.0500000 1026.47
1: 405050.97: 341.603 10.5350 0.0508866 27738.3 419.883 10.7618 0.0471600 2065.57 498.294 14.0557 0.0465954 1057.21
2: 115191.71: 341.507 10.5354 0.0556600 27858.3 421.299 10.1276 0.0418391 1986.87 503.484 13.9263 0.0356617 1262.92
3: 38076.077: 342.417 11.2347 0.0632863 27377.3 420.770 14.8188 0.0546385 2213.08 505.655 18.1187 0.0495791 1407.27
4: 36401.368: 343.360 11.7864 0.0723805 26889.9 426.228 23.2991 0.115937 2330.60 507.362 26.3221 0.0784007 1706.85
5: 16394.715: 343.437 11.7838 0.0741048 26883.4 423.172 19.5157 0.154983 2482.43 519.106 27.3302 0.165639 1558.34
6: 12437.878: 343.449 11.7884 0.0743107 26868.4 426.309 21.3703 0.207416 2569.34 517.635 24.8692 0.263019 1512.44
7: 12248.298: 343.429 11.7789 0.0740482 26885.0 426.114 20.9106 0.213771 2551.15 516.084 24.6528 0.200320 1527.81
8: 12235.845: 343.430 11.7791 0.0740580 26884.1 426.175 20.9728 0.214606 2555.89 515.690 24.4315 0.192340 1523.57
9: 12230.776: 343.430 11.7794 0.0740656 26883.7 426.227 20.9872 0.217407 2556.37 515.362 24.3697 0.180294 1523.84
10: 12217.446: 343.432 11.7803 0.0740936 26881.7 426.645 21.0955 0.238821 2558.55 514.148 24.1081 0.139162 1524.57
11: 12185.853: 343.435 11.7813 0.0741224 26879.7 427.203 21.2201 0.274725 2561.21 513.228 23.8124 0.126246 1525.05
12: 12174.404: 343.436 11.7819 0.0741410 26878.4 427.589 21.2985 0.310384 2564.07 512.065 23.4146 0.106315 1524.86
13: 12161.212: 343.437 11.7826 0.0741606 26877.1 427.933 21.3557 0.351018 2565.29 512.085 23.3748 0.109496 1524.09
14: 12155.955: 343.437 11.7826 0.0741621 26876.9 428.243 21.3974 0.394982 2565.96 511.729 23.2536 0.104486 1524.67
15: 12152.489: 343.438 11.7827 0.0741652 26876.7 428.497 21.4262 0.441353 2566.25 511.693 23.2270 0.104343 1524.34
16: 12150.125: 343.438 11.7829 0.0741713 26876.3 428.714 21.4500 0.491154 2566.61 511.637 23.2104 0.103651 1524.53
17: 12148.632: 343.438 11.7829 0.0741714 26876.3 429.008 21.4756 0.569129 2566.55 511.659 23.2185 0.103983 1524.51
18: 12147.015: 343.438 11.7827 0.0741674 26876.5 429.225 21.4869 0.653321 2566.19 511.648 23.2186 0.103855 1524.68
19: 12145.989: 343.438 11.7828 0.0741677 26876.4 429.391 21.4974 0.738613 2566.22 511.659 23.2218 0.103998 1524.65
20: 12145.369: 343.438 11.7829 0.0741710 26876.2 429.533 21.5074 0.830413 2566.43 511.651 23.2199 0.103902 1524.69
21: 12145.021: 343.438 11.7829 0.0741707 26876.2 429.685 21.5152 0.947698 2566.43 511.656 23.2211 0.103965 1524.66
22: 12144.714: 343.438 11.7828 0.0741698 26876.3 429.815 21.5202 1.08360 2566.35 511.653 23.2208 0.103927 1524.70
23: 12144.463: 343.438 11.7828 0.0741698 26876.3 429.913 21.5239 1.22124 2566.36 511.656 23.2217 0.103966 1524.69
24: 12144.317: 343.438 11.7829 0.0741705 26876.2 429.992 21.5273 1.35908 2566.42 511.651 23.2198 0.103907 1524.69
25: 12144.214: 343.438 11.7829 0.0741712 26876.2 430.059 21.5299 1.50140 2566.47 511.654 23.2205 0.103943 1524.67
26: 12144.204: 343.438 11.7829 0.0741712 26876.2 430.059 21.5300 1.51704 2566.50 511.650 23.2189 0.103890 1524.67
27: 12144.204: 343.438 11.7829 0.0741713 26876.2 430.059 21.5303 1.51705 2566.51 511.650 23.2188 0.103891 1524.67
28: 12144.204: 343.438 11.7829 0.0741714 26876.2 430.059 21.5305 1.51706 2566.53 511.651 23.2185 0.103891 1524.65
29: 12144.204: 343.438 11.7829 0.0741714 26876.2 430.059 21.5305 1.51706 2566.53 511.651 23.2185 0.103891 1524.65
30: 12144.204: 343.438 11.7829 0.0741714 26876.2 430.059 21.5305 1.51706 2566.53 511.651 23.2185 0.103891 1524.65
31: 12144.204: 343.438 11.7829 0.0741714 26876.2 430.059 21.5305 1.51706 2566.53 511.651 23.2185 0.103891 1524.65
I was having a similar problem so what I did was "force" new iterations by using the warnOnly=T
, which would result in actual estimates, then use those estimates as new starting values in a second nls()
. This is what my code ended up looking like:
a_start2 = 40
b_start2 = 200
p_start2 = 16.5 * mean(no.stage[Position == i & Stage == j & Year == k]) + 74.167
subset1 = which(Position == i & Stage == j & Year == k)
m2 = nls(
Percent ~ ((a) / sqrt(2*b*pi)) * exp(-(((DAFB - p)^2) / (2*b))),
start = list(a = a_start2, b = b_start2, p = p_start2),
control = list(maxiter = 50000, minFactor = 1/2000, warnOnly = TRUE),
algorithm = "port",
lower = list(a = 0.1, b = 100, p = -100),
upper = list(a = 200, b = 800, p = 400),
subset = subset1
)
print(summary(m2))
a_start3 = coef(summary(m2))["a", "Estimate"]
b_start3 = coef(summary(m2))["b", "Estimate"]
p_start3 = coef(summary(m2))["p", "Estimate"]
m3 = nls(
Percent ~ ((a) / sqrt(2*b*pi)) * exp(-(((DAFB - p)^2) / (2*b))),
start = list(a = a_start3, b = b_start3, p = p_start3),
control = list(maxiter = 50000, minFactor = 1/2000, warnOnly = TRUE),
algorithm = "port",
lower = list(a = 0.1, b = 100, p = -100),
upper = list(a = 200, b = 800, p = 400),
subset = subset1
)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With