0

この式に合わせる必要があります

 y ~ 1/(pi*a*(1+((x-2.15646)/a)^2))+1/(pi*b*(1+((x-2.16355)/b)^2))

変数 x と y に

 x<- c(2.15011, 2.15035, 2.15060, 2.15084, 2.15109, 2.15133, 2.15157, 2.15182, 2.15206, 2.15231, 2.15255, 2.15280, 2.15304, 2.15329, 2.15353, 2.15377, 2.15402, 2.15426, 2.15451, 2.15475, 2.15500, 2.15524, 2.15549, 2.15573, 2.15597, 2.15622, 2.15646, 2.15671, 2.15695, 2.15720, 2.15744, 2.15769, 2.15793, 2.15817, 2.15842, 2.15866, 2.15891, 2.15915, 2.15940, 2.15964, 2.15989, 2.16013, 2.16037, 2.16062, 2.16086, 2.16111, 2.16135, 2.16160, 2.16184, 2.16209, 2.16233, 2.16257, 2.16282, 2.16306, 2.16331, 2.16355, 2.16380, 2.16404, 2.16429, 2.16453, 2.16477, 2.16502, 2.16526, 2.16551, 2.16575, 2.16600, 2.16624, 2.16649, 2.16673, 2.16697, 2.16722, 2.16746, 2.16771, 2.16795, 2.16820, 2.16844, 2.16869, 2.16893, 2.16917, 2.16942, 2.16966, 2.16991)
 y<- c(3.77212,  3.79541,  3.84574,  3.91918, 4.01056,  4.11677,  4.23851,  4.37986, 4.54638,  4.74367,  4.97765,  5.25593,  5.58823,  5.98405,  6.44850,  6.98006, 7.57280,  8.22085,  8.92094,  9.66990, 10.45900, 11.26720, 12.05540, 12.76920, 13.34830, 13.74250, 13.92420, 13.89250, 13.67090, 13.29980, 12.82780, 12.30370, 11.76950, 11.25890, 10.80020, 10.41860, 10.13840,  9.98005,  9.95758, 10.07690, 10.33680, 10.73210, 11.25730, 11.90670, 12.67240, 13.54110, 14.49530, 15.51670, 16.58660, 17.67900, 18.75190, 19.74600, 20.59680, 21.24910, 21.66800, 21.83910, 21.76560, 21.46020, 20.94020, 20.22730, 19.35360, 18.36460, 17.31730, 16.26920, 15.26920, 14.35320, 13.54360, 12.85230, 12.28520, 11.84690, 11.54040, 11.36610, 11.32130, 11.39980, 11.59230, 11.88310, 12.25040, 12.66660, 13.09810, 13.50220, 13.82580, 14.01250)

x と y に従って 'a' と 'b' の値を推定します。'a' と 'b' は 0 から 1 の範囲です。

ただし、nls コマンドを使用すると、次のようになります。

nls(y ~1/(pi*a*(1+((x-2.15646)/a)^2))+1/(pi*b*(1+((x-2.16355)/b)^2)), control = list(maxiter = 500), start=list(a=0.4,b=0.4))

コンソールは次のエラーを報告しました:

singular gradient

コンソールにこのメッセージが表示される理由を誰か説明してもらえますか?

4

2 に答える 2

6

これにより、より適切な適合が得られます。

コード (以下) に入る前に、モデルにはいくつかの問題があります。

  1. これがプロトン NMR であると仮定すると、ピークの下の領域はプロトンの存在量 (つまり、プロトンの数) に比例します。モデルではこれが許可されておらず、基本的にすべてのピークが同じ面積を持つように強制されています。これがフィット感が悪い主な理由です。各ピークに「高さ」係数を含めることで、これに簡単に対応できます。
  2. モデルはピーク位置を想定しています。アルゴリズムに真のピーク位置を見つけさせないのはなぜですか?
  3. モデルはベースライン ドリフトを考慮していません。これは、ご覧のとおり、データセットでは非常に深刻です。モデルに線形ドリフト関数を追加することで、これに対応できます。

nls(...)このタイプのモデリングには適していません。使用するアルゴリズムは特に堅牢ではありません。デフォルトのアルゴリズムである Gauss-Newton は、オフセット データをフィッティングする場合に特に不十分です。そのため、ほぼ常に失敗するモデルで推定p1します。p2f(x-p1,x-p2)

nls.lm(...)より良いアプローチは、パッケージに実装されている非常に堅牢なレーベンバーグ・マルカート アルゴリズムを使用することminpackです。このパッケージは少し使いにくいですが、 ではアクセスできない問題を処理できますnls(...)。これを何度も行う場合は、ドキュメントを読んで、この例がどのように機能するかを理解する必要があります。

最後にnls.lm(...)、出発点も合理的でなければなりません。あなたのモデルでは a と b はピーク幅です。明らかに、それらはピーク位置の差に匹敵するか、それよりも小さくなければなりません。さもないと、ピークがくっつきます。(a,b) = (0.4, 0.4) の見積もりは大きすぎました。

plot(x,y)
library(minpack.lm)
lorentzian <- function(par,x){ 
  a  <- par[1]
  b  <- par[2]
  p1 <- par[3]
  p2 <- par[4]
  h1 <- par[5]
  h2 <- par[6]
  drift.a <- par[7]
  drift.b <- par[8]
  h1/(pi*a*(1+((x-p1)/a)^2))+h2/(pi*b*(1+((x-p2)/b)^2)) + drift.a + drift.b*x
}

resid   <- function(par,obs,xx) {obs-lorentzian(par,xx)}
par=c(a=0.001,b=0.001, p1=2.157, p2=2.163, h1=1, h2=1, drift.a=0, drift.b=0)
lower=c(a=0,b=0,p1=0,p2=0,h1=0, h2=0,drift.a=NA,drift.b=NA)
nls.out <- nls.lm(par=par, lower=lower, fn=resid, obs=y, xx=x, 
                  control = nls.lm.control(maxiter=500))
coef(nls.out)
#             a             b            p1            p2            h1            h2       drift.a       drift.b 
#  1.679632e-03  1.879690e-03  2.156308e+00  2.163500e+00  4.318793e-02  8.199394e-02 -9.273083e+02  4.323897e+02 

lines(x,lorentzian(coef(nls.out), x), col=2, lwd=2)

最後にもう 1 つ: SO の慣例では、回答を「受け入れる」まで 1 日待ちます。その理由は、受け入れられた回答を持つ質問が追加の注目を集めることはめったにないためです。回答を受け入れると、他の誰もそれを見ることはありません.

于 2014-01-29T19:54:51.490 に答える