RGRからロジスティック曲線を推定する

植物の重量に上限があったとして、相対成長速度が上限と今の重量の差に比例するとしてみよう。
${\displaystyle r \propto w_{max} - w }$
これは適当な定数をもってきて書き換えると、要するに相対成長速度が重量の一次式で表現できるということになる。
${\displaystyle r = a - bw }$
これを相対成長速度の定義に代入する。
${\displaystyle \frac{dw}{dt} = (a - bw)w }$
そして $w$ について解くとこうなる。
${\displaystyle w_t = \frac{a \exp (ac_1+a_t)}{b \exp (ac_1 + at) + 1} }$
これは実はロジスティック方程式である(定数を整理してみると良い)。

$a$ と $b$ と $c_1$ に適当な値を入れてプロットしてみよう。

a <- 0.4
b <- 0.2
c1 <- -15
t <- seq(0, 50, 1)

d <- logistic(t, c1, a, b)
plot(t, d)

f:id:Rion778:20161022203447p:plain
そして、重量と相対成長速度の関係は当然だが一次式になっている。

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(d))
r <- NA
for(i in 1:length(d)){
  r[i] <- r_t(d[c(i-1, i, i+1)], t[c(i-1, i, i+1)])
}
plot(d, r)

f:id:Rion778:20161022203615p:plain
点が不等間隔に並んでいるのも分かるだろう。重量の上限値でもって相対成長速度は0になるので、負の値になることはない。

ここからパラメータ $a$ と $b$ を求めるのは簡単で、単に回帰すれば良い。また $c_1$ はoptimizeを使えば求められる。ただし、 $c_1$ は曲線の立ち上がりを左右に移動させるだけなので、データのどこを $t=0$ と定めるかによって左右される。本質的なのは $a$ と $b$ だ。

plot(t, d)
result <- lm(r ~ d)
result_coef <- coef(result)
result_coef # > 0.4074815  -0.2040751 
f <- function(c1) (sum(logistic(t, c1, result_coef[1], -result_coef[2]) - d))^2
c1 <- optimize(f, interval = c(-100, 100))$minimum
points(t, logistic(t, c1, result_coef[1], -result_coef[2]), col = "red", type = "l")

f:id:Rion778:20161022204522p:plain

ゴンペルツ曲線の例との違い

今回の例では相対成長速度と時間を切り離している。そのため、例えば測定データのセットが複数あった場合、相対成長速度vs重量というデータに変換してしまうことで複数の結果を束ねて解析することが容易になる。

複数の葉や果実の生長、肥大を測定した場合に、どこを $t=0$ として測定を開始するのかというのは難しい問題である。似たような大きさのものを選んだとしても、特にサイズが小さいうちは指数関数的に生育が進むので、初期の誤差が生育曲線の位相の大きなずれとなって現れてくる。この問題を棚上げしたまま共通のパラメータを推定できるのは大きなメリットだろう。

ちなみにゴンペルツ曲線とロジスティック曲線はそもそも形が違う(ゴンペルツ曲線は対称ではない)ので、その点は注意が必要である。