曲線のフィッティング

データがある曲線 y = a * ( 1 – exp( -b * t ) ) に従っているとします。
a と b は関数のパラメータです。
例えば a = 2, b = 0.7 の場合には、この曲線は次のようになります。

library( ggplot2 )
library( ggthemes )

# 曲線の関数
mycurve = function( x, a, b ){
  return( a * ( 1 - exp( -b * x ) ) )
}

# データの作成
a = 2
b = .7
mydata = data.frame(
  x = c ( 0.1, 0.2, 0.5, 1, 3, 6 ),
  y = mycurve( c( 0.1, 0.2, 0.5, 1, 3, 6 ), a, b )
)

p = ggplot( mydata, aes( x = x, y = y ) ) +
  stat_function( fun = mycurve, args = list( a = a, b = b ), size = 1.5 ) +
  geom_point( size = 4 ) +
  theme_bw()

plot( p )

optim関数によるパラメータ推定

次に、この関数に従っていると想定されるデータがある時に、関数のパラメータを推定します。
推定には optim 関数を使用します。
optim 関数には推定したいパラメータの初期値と、最小化したい関数を指定します。
最小二乗法を使う場合は、最小化したい関数は二乗誤差になります。
以下の例では最小化したい関数を errorfn という名前の関数で作成しています。

library( ggplot2 )
library( ggthemes )

# 曲線の関数
mycurve = function( x, a, b ){
  return( a * ( 1 - exp( -b * x ) ) )
}

# データの作成
a = 2
b = .7
mydata = data.frame(
  x = c (0.1, 0.2, 0.5, 1, 3, 6 ),
  y = mycurve( c(0.1, 0.2, 0.5, 1, 3, 6 ), a, b )
)

# データにノイズを加える
mydata$y = mydata$y + rnorm( 6, mean = 0, sd = .1 )

# 最小二乗法によるパラメータの推定( p はパラメータ、d はデータです)
errorfn = function( p, d ) {
  y = sum( ( mycurve( d$x, p[ 1 ], p[ 2 ]) - d$y )^2 )
}
params = optim( par = c( 1, 1 ), fn = errorfn, d = mydata )$par
# par は推定したいパラメータの初期値で、適当な値を指定します(ここでは a と b の両方とも1)
# fn は最小化したい関数で、この関数の引数 d に mydata を指定しています

print( params )

# 結果の図示
p = ggplot( mydata, aes( x = x, y = y ) ) +
  stat_function( fun = mycurve, args = list( a = a, b = b ), size = 1.5 ) +
  stat_function( fun = mycurve, args = list( a = params[1], b = params[2] ), color = "red", size = 1.5 ) +
  geom_point( size = 4 ) +
  theme_bw()

plot( p )

params という変数にパラメータの推定結果が入っています。
a = 2.0209496, b = 0.795333 という結果が得られました。
a = 2, b = 0.7 が正解なので、近い値が出ています。
(データにノイズを加えずにやればほぼ正解と同じ値になります)
赤線は推定されたパラメータを使って引いた(データにフィットさせた)曲線です。
黒線はノイズを含まない場合の曲線(真のモデル)です。

コメント