Excel などで（もちろん R ででも），初期値をいろいろ変えてあてはめ曲線がどのように変化するか，試行錯誤するのも基本的に必要でしょう。しかし，ある程度もっともらしい初期値が見つかったら，さらによい初期値を探すには，optim という関数があります。具体例を挙げて説明しましょう。

例： 表 1 のような x，y の測定値に対して，y = α / { 1 + β exp(-γx) } + δ という関数をあてはめ，パラメータ，α，β，γ，δを求める。

表 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)

# あてはめるべき関数 a/(1+b*exp(-c*x))+d で予測値を計算し，実測値との差の二乗和（残差平方和）を計算する関数を定義する
resid <- function(par)	# 関数名は何でも良い。引数は，パラメータのベクトル。x, y は，大域変数として参照する。
{
	yhat <- par[1]/(1+par[2]*exp(-par[3]*x))+par[4]	# y ~ a/(1+b*exp(-c*x))+d のこと
	sum((y-yhat)^2)	# 残差平方和を返す
}

optim の起動と結果の出力	# パラメータの詳細は ? optim で

> optim(c(-100, 50, 1, 100), resid)	# 適当な初期値ベクトルを第1引数に指定
$par
[1] -82.29288172  42.99372287   0.09201423 100.99425999	# 解

$value
[1] 18.91874	# そのときの関数値（残差平方和）

$counts
function gradient 
     501       NA		# 関数が呼ばれた回数

$convergence
[1] 1	# 0 のときは収束した。1 ならまだ収束しなかった。

$message
NULL

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