読者です 読者をやめる 読者になる 読者になる

RGRからゴンペルツ曲線を推定する 3

R

ゴンペルツ曲線では相対成長速度と重量の関係はどうなってるのかと思って適当に回帰してみたら

r = a - b\ln w

が良く当てはまったので、この仮定(相対成長速度は重量の対数に比例して減少する)のもとに微分方程式

\frac{1}{w}\frac{dw}{dt} = a - b \ln w

を解いてみると、ちゃんとゴンペルツ関数が出て来る。

w_t = c_1c_2^{\exp (-bt)}

ここで

c_1 = \exp(a/b)

c_2 = \exp(\exp(bc)/b)

であり、cに時間に関するパラメータが分離されているので、始点の問題を棚上げできる。

あとは適当に。

gompertz <- function(t, c1, a, b){ # RGR vs Time
  c1 * exp(-(exp(a-b*t)/b))
}
t <- seq(0, 20, 0.5)
a <- 0.3
b <- 0.25
c1 <- 10
g <- gompertz(t, c1, a, b)
plot(t, g)

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)])
}

gompertz2 <- function(t, c1, a, b){ # RGR vs Weight
  exp(a/b - exp(b*c1 - b*t)/b)
}
result <- lm(r ~ log(g))
result_coef <- coef(result)
f <- function(c1) (sum(gompertz2(t, c1, result_coef[1], -result_coef[2]) - g))^2
c1 <- optimize(f, interval = c(-50, 100))$minimum
points(t, gompertz2(t, c1, result_coef[1], -result_coef[2]), col = "red", type = "l")

f:id:Rion778:20161023090414p:plain