1

を使用してプロットに最適な行を追加する基本的な関数を作成しようとしていますnls。に渡された数式によってデータが正確に定義されていない限り、これはうまく機能しnlsます。私は問題を認識しており、これはここで報告されているように文書化された動作であることを認識しています。

私の質問は、モデルによって正確に記述されているデータに関係なく、どうすればこれを回避し、最適な線を強制的にプロットできるかということです. データの一致を正確に検出し、完全に適合する曲線をプロットする方法はありますか? 私の現在の危険な解決策は次のとおりです。

#test data
x <- 1:10
y <- x^2
plot(x, y, pch=20)

# polynomial line of best fit
f <- function(x,a,b,d) {(a*x^2) + (b*x) + d}
fit <- nls(y ~ f(x,a,b,d), start = c(a=1, b=1, d=1)) 
co <- coef(fit)
curve(f(x, a=co[1], b=co[2], d=co[3]), add = TRUE, col="red", lwd=2)

次のエラーで失敗します。

Error in nls(y ~ f(x, a, b, d), start = c(a = 1, b = 1, d = 1)) : 
  singular gradient

私が適用する簡単な修正jitterはデータにわずかですが、これは少し破壊的でハックのようです。

# the above code works after doing...
y <- jitter(x^2)

より良い方法はありますか？

4

1 に答える 1

6

Levenberg-Marquardt を使用します。

x <- 1:10
y <- x^2

f <- function(x,a,b,d) {(a*x^2) + (b*x) + d}
fit <- nls(y ~ f(x,a,b,d), start = c(a=1, b=0, d=0)) 

Error in nls(y ~ f(x, a, b, d), start = c(a = 1, b = 0, d = 0)) : 
  number of iterations exceeded maximum of 50

library(minpack.lm)
fit <- nlsLM(y ~ f(x,a,b,d), start = c(a=1, b=0, d=0))
summary(fit)

Formula: y ~ f(x, a, b, d)

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a        1          0     Inf   <2e-16 ***
b        0          0      NA       NA    
d        0          0      NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0 on 7 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 1.49e-08

開始値を調整する必要があり、結果は開始値に敏感であることに注意してください。

fit <- nlsLM(y ~ f(x,a,b,d), start = c(a=1, b=0.1, d=0.1))

Parameters:
    Estimate Std. Error    t value Pr(>|t|)    
a  1.000e+00  2.083e-09  4.800e+08  < 2e-16 ***
b -7.693e-08  1.491e-08 -5.160e+00  0.00131 ** 
d  1.450e-07  1.412e-08  1.027e+01  1.8e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 6.191e-08 on 7 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 1.49e-08
于 2012-12-19T08:13:50.593 に答える