8

このトピックに関する他のスレッドを検索してみましたが、どの修正も機能していません。自然実験の結果があり、イベントの連続発生回数が指数分布に適合することを示したいと考えています。私のRシェルは以下に貼り付けられています

f <- function(x,a,b) {a * exp(b * x)}
> x
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
[26] 26 27
> y
 [1] 1880  813  376  161  100   61   31    9    8    2    7    4    3    2    0
[16]    1    0    0    0    0    0    1    0    0    0    0    1
> dat2
    x    y
1   1 1880
2   2  813
3   3  376
4   4  161
5   5  100
6   6   61
7   7   31
8   8    9
9   9    8
10 10    2
11 11    7
12 12    4
13 13    3
14 14    2
> fm <- nls(y ~ f(x,a,b), data = dat2, start = c(a=1, b=1)) 
Error in numericDeriv(form[[3L]], names(ind), env) : 
  Missing value or an infinity produced when evaluating the model
> fm <- nls(y ~ f(x,a,b), data = dat2, start = c(a=7, b=-.5)) 
Error in nls(y ~ f(x, a, b), data = dat2, start = c(a = 7, b = -0.5)) : 
  singular gradient
> fm <- nls(y ~ f(x,a,b), data = dat2, start = c(a=7,b=-.5),control=nls.control(maxiter=1000,warnOnly=TRUE,minFactor=1e-5,tol=1e-10),trace=TRUE) 
4355798 :   7.0 -0.5
Warning message:
In nls(y ~ f(x, a, b), data = dat2, start = c(a = 7, b = -0.5),  :
  singular gradient

書式設定が悪いことをお許しください。最初の投稿はこちらです。x にはヒストグラムのビンが含まれ、y にはそのヒストグラム内の各ビンの出現回数が含まれます。dat2 は 14 で途切れます。なぜなら、0 カウントのビンは指数関数的回帰を放棄するからです。実際には最初の 14 だけを適合させる必要があります。14 を超えるカウントを持つビンは、それらが特別であると信じる生物学的理由があります。私が最初に得た問題は無限大でしたが、値が 0 ではないため得られませんでした。ここの別の投稿で示唆されているように適切な開始値を指定すると、特異な勾配エラーが発生します。私が見た唯一の他の投稿には、より多くの変数があり、反復回数を増やしてみましたが、成功しませんでした。どんな助けでも大歓迎です。あ

4

2 に答える 2

19

1) 線形化して開始値を取得する より良い開始値が必要です。

# starting values
fm0 <- nls(log(y) ~ log(f(x, a, b)), dat2, start = c(a = 1, b = 1))

nls(y ~ f(x, a, b), dat2, start = coef(fm0))

与える:

Nonlinear regression model
  model: y ~ f(x, a, b)
   data: x
        a         b 
4214.4228   -0.8106 
 residual sum-of-squares: 2388

Number of iterations to convergence: 6 
Achieved convergence tolerance: 3.363e-06

1a)同様に、次lmのように記述して初期値を取得するために使用できます。

y ~ a * exp(b * x)

なので

y ~ exp(log(a) + b * x)

両方のログを取得して、log(a) と b で線形のモデルを取得します。

log(y) ~ log(a) + b * x

これは次を使用して解決できますlm。

fm_lm <- lm(log(y) ~ x, dat2)
st <- list(a = exp(coef(fm_lm)[1]), b = coef(fm_lm)[2])
nls(y ~ f(x, a, b), dat2, start = st)

与える:

Nonlinear regression model
  model: y ~ f(x, a, b)
   data: dat2
       a        b 
4214.423   -0.811 
 residual sum-of-squares: 2388

Number of iterations to convergence: 6 
Achieved convergence tolerance: 3.36e-06

1b)再パラメータ化することで動作させることもできます。その場合、パラメータ変換に合わせて初期値を変換すれば、a = 1 と b = 1 が機能します。

nls(y ~ exp(loga + b * x), dat2, start = list(loga = log(1), b = 1))

与える:

Nonlinear regression model
  model: y ~ exp(loga + b * x)
   data: dat2
  loga      b 
 8.346 -0.811 
 residual sum-of-squares: 2388

Number of iterations to convergence: 20 
Achieved convergence tolerance: 3.82e-07

したがって、b は示されているとおりであり、a = exp(loga) = exp(8.346) = 4213.3

2) plinear より簡単なもう 1 つの可能性は、使用することalg="plinear"です。この場合、線形に入るパラメーターに開始値は必要ありません。その場合、b=1質問の開始値で十分と思われます。

nls(y ~ exp(b * x), dat2, start = c(b = 1), alg = "plinear")

与える:

Nonlinear regression model
  model: y ~ exp(b * x)
   data: dat2
        b      .lin 
  -0.8106 4214.4234 
 residual sum-of-squares: 2388

Number of iterations to convergence: 11 
Achieved convergence tolerance: 2.153e-06
于 2013-08-21T18:06:22.660 に答える