Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

nls() false convergence (despite good starting values)

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
like image 325
Ryan Avatar asked Jun 27 '16 23:06

Ryan


1 Answers

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
)
like image 194
Al Kovaleski Avatar answered Sep 29 '22 15:09

Al Kovaleski