前回はRGRを重量の対数の変化速度として求めたが、そもそもRGRの定義は
であったので、複数のデータポイントがあれば何らかの方法で平滑化して微分すれば、より精度良くRGRが推定できるだろう。例えば前後2点を含めた3点をラグランジュ補間して得た二次関数を微分するとか。
Rでラグランジュ補間をするパッケージは多分あるが面倒なので二次関数を当てはめてしまう。
r_t <- function(w, t) sum(coef(lm(w ~ t + I(t^2)))[2:3] * c(1, 2*t[2])) / w[2] r <- numeric(length(g)) r <- NA for(i in 1:length(g)){ r[i] <- r_t(g[c(i-1, i, i+1)], t[c(i-1, i, i+1)]) }
そしてあとは前回のヤツを多少いじって実行する。
result <- glm(r ~ t, family = gaussian(link="log"), start = c(1, 0)) result_coef <- coef(result) plot(t, g) f <- function(c1) (sum(gompertz(t, c1, result_coef[1], -result_coef[2]) - g))^2 c1 <- optimize(f, interval = c(0, 100))$minimum points(t, gompertz(t, c1, result_coef[1], -result_coef[2]), col = "red", type = "l")
すると少しはマシな推定値が出る。
ただ、繰り返しになるがゴンペルツ曲線を求めるなら最初からゴンペルツ曲線をデータに当てはめたほうが良い。例えばgに適当な誤差を与えて推定を繰り返してみると、マズさが分かるだろう。
特定のタイミングにおけるRGRの推定値がほしければ、ゴンペルツ関数を微分したものに推定したパラメータを渡せば良い。微分した結果を関数として使うにはderiv()を使って引数にfunc=TRUEを指定する。
dg <- deriv(~ c1 * exp(-(exp(a-b*t)/b)), "t", func = TRUE) attr(dg(t), "gradient") r <- attr(dg(t), "gradient")[,1] / gompertz(t, c1, a, b) plot(r)
