2016-06-27 9 views
4

私は3つの指数関数的に修正されたガウス分布(EMG)を畳み込み曲線にフィットさせる曲線適合スクリプトに取り組んでいます。私の基本関数はガウス分布に似ていますが、関数の指数関数の成分の重みを決定する第3のパラメータ(最初の2つはmusigmaです)を含みます。nls()false収束(良い開始値にもかかわらず)

したがって全体、各EMGピークはEMGが最小化されるパラメータの総数をピーク3をデコンボリューションする

(値> 1.0と実験データと一致するために)3つのパラメータ、プラス振幅係数を取り3x4のは= 12

いくつかのケースではフィットがうまく動作しますが、多くの場合、それが収束に失敗し、わずか50かそこらの反復(​​井戸bの後、この

Convergence failure: false convergence (8) 

のようなエラーを返します限度額を超えてください)。

トレースオプションを使用すると、結果がデータの照合に非常に近いことがわかります。

データ=黒(ノイズを添加)、初期=オレンジ、最終:と私の初期推定値から曲線をプロットすることによって、また、起動用パラメータがデータに合理的な近接範囲内にあることがわかります

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", 
) 

アルゴリズムが成功しているように見えるときに、なぜ私はエラーを取得しています:エラーここで=赤

前の反復は、私がnls()を呼ぶ私のコードのサンプルを、ありますか?

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 
+0

(:?これは私のコードは次のように見てしまったものです「コントロール引数でwarnOnlyは= TRUEの設定Rバージョン2.5以降(非集中型のオブジェクトを返します(nls.controlを参照してください) .0)これはさらなる収束解析には有用かもしれないが、推論のためには有用ではないかもしれない。 –

+0

パラメータに応じてSSEをプロットします。等高線プロットを使用すると、2つのパラメータの変動に対するSSEの感度を調べることができます。私は、明確な最小値がないことが分かると思います。 – Roland

+0

ええと残念ながら私は12のパラメータを持っているので、SSEとの関係を視覚化する可能性はほとんどありません。しかしSSEが最小に達しているので、私はあなたの意見を見ますが、アルゴリズムはSSEの結果変化なしでパラメータを修正し続け、誤った収束を宣言します。ステップサイズが小さすぎるかもしれませんか?トレース出力の私の更新されたポストを参照してください – Ryan

答えて

2

/私はこれをコメントしようとしましたが、私はちょうど私のプロファイルを作成しました。

とにかく私は、warnOnly=Tを使って新しい反復を強制して、実際の見積もりを得て、その推定値を2番目の新しい開始値として使用するようにしました。nls() `)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) 
関連する問題