nls という関数があります。具体例を挙げて説明しましょう。

例： 表 1 のような x，y の測定値に対して，y = α / { 1 + β exp(-γx) } + δ という関数をあてはめ，パラメータ，α，β，γ，δを求める（同じ問題を Excel で解く方法）。

表 1. 二変数データ
   x    y
   0   98
   5   97
  10   96
  15   94
  20   92
  25   84
  30   78
  35   73
  40   61
  45   53
  50   43
  55   37
  60   30
  65   28
  70   24
  75   22
  80   21
  85   20
  90   20
  95   20
 100   18

# x は 0 から始まり，5 刻みで，100 までの 21 個の値を取る
> x <- seq(0,100,5)
# y は，順に 21 個の値を取る
> y <- c(98, 97, 96, 94, 92, 84, 78, 73, 61, 53, 43, 37, 30, 28, 24, 22, 21, 20, 20, 20, 18)
# データフレーム dat を作る
> dat <- data.frame(x, y)
> dat # 作ったのは，以下のようなデータ
     x    y
1    0   98
2    5   97
3   10   96
： 【中略】
19  90   20
20  95   20
21 100   18
# y ~ a/(1+b*exp(-c*x))+d が，あてはめるべき関数
# a, b, c, d はパラメータ（これを推定する）
# start= で，初期値を指定する 初期値の選び方
> nls(y ~ a/(1+b*exp(-c*x))+d, dat, start=list(a=1, b=1, c=1, d=1))
# 初期値が不適切な場合には，以下のようなエラーメッセージが出るので，初期値を変えてみる
Error in nls(y ~ a/(1 + b * exp(-c * x)) + d, dat, start = list(a = 1,  :
	singular gradient
# nls 関数の引用だけだと，モデル関数，推定されたパラメータ値，残差平方和が出力される
> nls(y ~ a/(1+b*exp(-c*x))+d, dat, start=list(a=-80, b=47, c=0.1, d=100))
Nonlinear regression model
  model:  y ~ a/(1 + b * exp(-c * x)) + d
   data:  dat
           a            b            c            d
-81.19250400  47.73105651   0.09353186 100.00271945
 residual sum-of-squares:  16.63448
# nls 関数の結果をオブジェクトに付値（代入）すると，後で，いろいろな結果を選択して出力することができる
> result <- nls(y ~ a/(1+b*exp(-c*x))+d, dat, start=list(a=-100, b=50, c=0.2, d=100))
# オブジェクトに対して sumamry 関数を適用すると，パラメータの推定値，標準誤差，t 値，P 値，のほかに，残差標準誤差，パラメータ間の相関係数が出力される
> summary(result)

Formula: y ~ a/(1 + b * exp(-c * x)) + d

Parameters:
   Estimate Std. Error t value Pr(>|t|)   
a -81.19250    1.03708 -78.290  < 2e-16 ***
b  47.73111    6.52126   7.319 1.20e-06 ***
c   0.09353    0.00300  31.180  < 2e-16 ***
d 100.00272    0.76138 131.343  < 2e-16 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

Residual standard error: 0.9892 on 17 degrees of freedom

Correlation of Parameter Estimates:
        a       b       c
b  0.8328               
c  0.8249  0.9709       
d -0.8968 -0.8118 -0.7241

直前のページへ戻る　　 E-mail to Shigenobu AOKI