S字?型の形状を取る。 ゆっくり立ち上がり、ググっと上がって、飽和してく。 耐久消費財の普及率が本では例示されてるけど、特定集団内での何らかの変化を描写するのに向いてるはず。
dy/dt = r * (K - y) * y
定義?としては、こういう形。
微分を解いて、 (Y / (K - Y) ) = C * exp(krt)
Yについて解いて Y = (K / (1 - (1/C) * exp(-rkt) )
tが大きくなれば、分母は1になってく。なので、YはKに近づいていく。
また、d^2(x)/ (dy)^2 二次微分をとると、 d^2 (x)/ (dy)^2 = r^2 * K * (K-Y) * (K - 2Y)
Y = K/2 のところにくると、変曲点を迎える。 勢いが変化し、終盤戦に入るイメージか?
下で行くと、そのまんまなのは、TVくらいかな。(Rで学ぶ時系列入門のデータを打ち込んだもの)
形状式の部分の、(1/C) expの係数と、expの中を変えて、形状を見る。 R> curve(100/(1+100*exp(-1*0.2*x)),xlim=c(1,100), main="basic")
R> curve(100/(1+100*exp(-1*0.5*x)),xlim=c(1,100), main="増加早い")
R> curve(100/(1+500*exp(-1*0.2*x)),xlim=c(1,100), main="立ち上がりゆっくり")
R> curve(100/(1+500*exp(-1*0.5*x)),xlim="両方足す")
推定(あてはめ)一次式に変換したり、変化率から推定したりがあるみたいだけど、飛ばして、 Rの非線形最小二乗法の関数を使う。 nls
Rで学ぶ時系列入門では、微分式が一次式なので、そこから近似してるけど、スキップする。 この二つを使っていく。
初期値の設定が問題。それっぽいのを見つけて、最適化する必要がある。 また、推定方法も誤差のニ乗和ではなく、 誤差式の偏微分がゼロ、極小値を求めてるだけなので、なおさら重要? メンドウだ、、、
とりあえず、今回のデータは、40年くらいで100%近くまで登っていくデータなので、 #最適化する関数を定義
> resid <- function(prm){
+ rhat <- prm[1]/(1+prm[2] * exp(-prm[3]*time)) + prm[4]
+ sum((rate - rhat)^2)
+ }
#初期値を最適化する関数を呼ぶ。初期値は、、、 opt <- optim( c(100,100,0.5,0), resid) R> opt $par [1] 103.5381 15.0417 0.1626 -0.5094 $value [1] 413.7 $counts function gradient 501 NA $convergence [1] 1 $message NUL #この初期値をもとに、非線形最小二乗法関数を呼ぶ R> rt <- nls(rate ~ a/(1+b*exp(-c*time))+d, start=list(a=opt$par[1],b=opt$par[2],c=opt$par[3],d=0)) 値を抽出 R> pr <- coef(rt) 実績値をplot R> plot(rate) 推計値をlinesで描写 R> lines(pr[1] / (1 + pr[2] * exp(-1 * pr[3] * time)) + pr[4])
当てはめの方は、ラインの伸ばせば、未来に対する推計になる。もっとも、上限値が(103 - 0.5)と、100%をオーバーしてるけど、、、パラメータで、上限は固定にするべきだった、、、 |