nls()错误收敛(尽管起始值很好)

2023年5月1日 353次阅读

我一直致力于曲线拟合脚本,它将3个指数修正的高斯(EMG)拟合成卷积曲线.我的基函数类似于高斯分布,但包括第三个参数(前两个是mu和sigma),它确定函数的指数分量的权重.

总的来说,每个EMG峰值需要3个参数,加上幅度系数(为了使实验数据与值> 1.0匹配)

为了解卷积3个EMG峰值,要最小化的参数总数为3×4 = 12

在某些情况下,拟合很好地工作,但在许多情况下它无法收敛,并返回这样的错误

Convergence failure: false convergence (8)

经过50次左右的迭代(远低于极限).

通过使用trace选项,我可以看到结果与数据匹配非常接近.通过绘制初始估算值的曲线,还可以看出起始参数与数据的合理接近程度：

数据=黑色(添加噪声),初始=橙色,错误前的最终迭代=红色

这是我的代码示例,我在其中调用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",
)

那么为什么我在算法似乎成功时遇到错误？

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

最佳答案 /我试图评论这个,但我刚创建了我的个人资料,所以我不能.

无论如何,我遇到了类似的问题所以我所做的是通过使用warnOnly = T“强制”新的迭代,这将导致实际估计,然后在第二个nls()中使用这些估计作为新的起始值.这就是我的代码最终看起来像：

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=T),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=T),algorithm ="port", lower=list(a=0.1, b=100, p=-100), upper=list(a=200, b=800, p=400), subset=subset1)
点赞