成長曲線をpythonで描いてみよう

# 目次
[趣旨](#趣旨)
[vonBentalanffy](#vonbentalanffy)
[Gompertz](#gompertz)
[Logistic](#logistic)
[Richards](#richards)
[追記](#追記)
[参考](#参考)

# 趣旨
講義にて、水産資源学における色々な成長曲線を学んだので、pythonの練習も兼ねて作画してみようと思います。
なお今回はpythonの使用が目的なので、詳細な水産資源学の内容は専門書等をご覧ください。

今回使うpythonのライブラリは``matplotlib.pyplot``と``numpy``です。
標準で装備されている``math``ライブラリも使用しています。

リストの内包表現を用いて計算し、その結果を図示します。
内包表記については以下のサイトなどを併せてご覧ください。

https://aiacademy.jp/media/?p=1252

この記事は一部、説明が丁寧でない箇所があります(プログラミングや指揮に関する詳細な記述がないです、サンプルコード程度にご高覧ください)。
気づいたら直していこうと思います()

## vonBentalanffy
von Bentalanffyの成長曲線は以下の式で表される式です。

$l_{t} = l_{\infty}(1-e^{-K(t-t0)})$

$l_{t}$：t歳の体長
$l_{\infty}$：極限体長
$K$：成長係数
$t_{0}$：原点を調整するためのパラメータ

以下、pythonにおける作画のコードです。

```python
t = list(range(0,25,1))
l_inf = 20 # 極限体長
K = 0.3 # 成長係数

lt = [l_inf*(1-math.e**(-K*i)) for i in t]

plt.plot(t,lt)
plt.show()
```

結果
![von_Bentalanffy.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3639393/be81de26-c592-d3a2-4475-384162aabd43.png)

では次に成長係数を変化させて図を描いてみます。
成長係数が0.3, 0.5, 0.7の3つの場合の曲線を比較してみます。

```python
t = list(range(0,25,1))
l_inf = 20

K1 = 0.3
K2 = 0.5
K3 = 0.7

lt1 = [l_inf*(1-math.e**(-K1*i)) for i in t]
lt2 = [l_inf*(1-math.e**(-K2*i)) for i in t]
lt3 = [l_inf*(1-math.e**(-K3*i)) for i in t]

fig, ax = plt.subplots()
ax.set_xlabel('t')  # x軸ラベル
ax.set_ylabel('lt')  # y軸ラベル
ax.grid() # 罫線を引く
ax.plot(t, lt1, color="black", label="K=0.3") # 作画
ax.plot(t, lt2, color="red", label="K=0.5") # 作画
ax.plot(t, lt3, color="blue", label="K=0.7") # 作画
ax.legend(loc=0)

plt.show()
```

結果
![von_Bentalanffy_multiple.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3639393/d29538ca-4b31-da0a-b989-a296626483f0.png)

# Gompertz
Gompertzの成長曲線は以下の式で表されます。

$l_{t} = l_{0}\exp(-e^{-K(t-t_{0})})$

$l_{t}$：t歳の体長
$l_{\infty}$：極限体長
$K$：成長係数
$t_{0}$：変曲点を表すパラメータ

pythonで作図してみましょう、適当な変数を代入して描いてみます。

```python
t = list(range(0,25,1))
l_inf = 20
K = 0.5
t0 = l_inf / math.e # 新たに追加、この式はGompertzの式を微分して導出できます

lt = [l_inf*math.e**(-math.e**(-K*(i-t0)) for i in t]

plt.plot(t,lt)
plt.show()
```

結果
![Gompertz.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3639393/907f8f22-539c-9ea8-975c-e91bc4cc730e.png)

Gompertz式に関して、y軸を対数軸に直してやるとvon Bentalanffy式の形になります。
気になる人は数式をいじってみてください。
以下ではコードと結果のみ示します。

```python
t = list(range(0,25,1))
l_inf = 20
K = 0.5
t0 = l_inf/math.e

lt = [l_inf*math.e**(-math.e**(-K*(i-t0))) for i in t]

fig, ax = plt.subplots()
ax.set_xlabel('t')
ax.set_ylabel('lt')
ax.set_yscale('log') # y軸を対数軸に変更
ax.grid()
ax.plot(t, lt)

plt.show()
```

結果
![Gompertz_longarithm.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3639393/c9bc11d6-ccda-cc78-c4e7-9f17681b689d.png)

Gompertz式に関しても複数の成長係数を比較してみましょう。
成長係数が0.3, 0.5, 0.7の3つの場合の曲線を比較してみます。

```python
t = list(range(0,25,1))
l_inf = 20
t0 = l_inf/math.e

K1 = 0.3
K2 = 0.5
K3 = 0.7

lt1 = [l_inf*math.e**(-math.e**(-K1*(i-t0))) for i in t]
lt2 = [l_inf*math.e**(-math.e**(-K2*(i-t0))) for i in t]
lt3 = [l_inf*math.e**(-math.e**(-K3*(i-t0))) for i in t]

fig, ax = plt.subplots()
ax.set_xlabel('t')
ax.set_ylabel('lt')
ax.grid()
ax.plot(t, lt1, color="black", label="K=0.3")
ax.plot(t, lt2, color="red", label="K=0.5")
ax.plot(t, lt3, color="blue", label="K=0.7")
ax.legend(loc=0)

plt.show()
```

結果
![Gompertz_multiple.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3639393/b40f00d0-3971-6da3-4a7e-626c1a6c943a.png)

7ぐらいのところに変曲点があることがわかります。
Gompertz式の変曲点は$\frac{l_\infty}{e}$となることが知られています(微分して導出可能)ので、良さげですね。

# Logistic

生態学では御用達、Logistic式です。
Logistic式は以下のように表されます。

$l_{t} = l_{\infty}\left(1+e^{-K(t-t_{0})} \right)^{-1}$

$l_{t}$：t歳の体長
$l_{\infty}$：極限体長
$K$：成長係数
$t_{0}$：変曲点を表すパラメータ

pythonで作図してみましょう、適当な変数を代入して描いてみます。

```python
l_inf = 20
K = 0.5
t_0 = l_inf/2

t = list(range(0,25,1))

lt = [l_inf * (1 + math.exp(-K * (ti - t_0))) ** (-1) for ti in t]

fig, ax = plt.subplots()
ax.set_xlabel('t')
ax.set_ylabel('lt')
ax.grid()
ax.plot(t, lt)

plt.show()
```

結果
![Logistic.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3639393/0bdc2ac3-cd7e-f0f2-8d7d-9a9e475c8e3e.png)

Gompertz式と見比べると...変曲点が少し右にずれていることに気づくかもしれません。
Logistic式の変曲点は$\frac{l_{\infty}}{2}$となることが数式から導出できます。

Logistic式についても複数の成長係数を比較してみましょう。
成長係数が0.3, 0.5, 0.7の3つの場合の曲線を比較してみます。

```python
K1 = 0.3
K2 = 0.5
K3 = 0.7

lt1 = [l_inf * (1 + math.exp(-K1 * (ti - t_0))) ** (-1) for ti in t]
lt2 = [l_inf * (1 + math.exp(-K2 * (ti - t_0))) ** (-1) for ti in t]
lt3 = [l_inf * (1 + math.exp(-K3 * (ti - t_0))) ** (-1) for ti in t]

fig, ax = plt.subplots()
ax.set_xlabel('t')
ax.set_ylabel('lt')
ax.grid()
ax.plot(t, lt1, color="black", label="K=0.3")
ax.plot(t, lt2, color="red", label="K=0.5")
ax.plot(t, lt3, color="blue", label="K=0.7")
ax.legend(loc=0)

plt.show()
```

結果
![Logistic_multiple.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3639393/deba62c4-62bc-2c11-e232-7b57c66aa704.png)

# Richards

Richards式は以下の式で表されます。

$l_{t} = l_{\infty}(1+re^{-K(t-t_{0})})^{-\frac{1}{r}}$

$l_{t}$：t歳の体長
$l_{\infty}$：極限体長
$K$：成長係数
$t_{0}$：変曲点を表すパラメータ

ここでrは任意の数です。
- r=-1の時に[von Bentalanffy式](#vonBentalanffy)
- r=0の時に[Gompertz式](#Gompertz)
- r=1の時に[Logistic式](#Logistic)
となります。

pythonで作画してみましょう。$l_{\infty}$と$K$は適当な値を代入しておきます。
$t_{0}$は0としておきます。

```python
t_value = np.arange(-2, 25, 0.1)  # 0から25の範囲を100点でサンプリング
t_value = list(np.round(t_value, 1))
t_value.remove(-0.)
l_inf = 20
K = 0.5
t_0 = 0

fig, ax = plt.subplots()
ax.set_xlabel('t')
ax.set_ylabel('lt')
ax.grid()

r_values = [-1, -1/3, 1, 2, 5, math.inf]  # 描画する r の値を指定

for r in r_values:
    lt_values = [l_inf * (1 + r * math.exp(-K * (t - t_0))) ** (-1/r) for t in t_value]
    ax.plot(t_value, lt_values, label=f"r={r}")

lt_values = [l_inf * math.exp(-math.e ** (-K * (t - t_0))) for t in t_value]
ax.plot(t_value, lt_values, label=f"r=0")

ax.legend(loc=0)
plt.show()
```

結果

![Richards.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3639393/30352f56-1070-9a31-ffb6-10e5563fc928.png)

確かにそれぞれの値で、それぞれの曲線に近しい形が出力されていることがわかります。

以上

# 追記

@Wolfmoon さんがnumpyを活用したコードの訂正を紹介してくださいました。

https://qiita.com/WolfMoon/items/6acbe1e2a6ff97e90dd3#_reference-4115a0eb727cc454288e

# 参考

https://aiacademy.jp/media/?p=1252

専門書の例

https://kaibundo.jp/2022/02/56334/

その他

http://katsuo.fs.a.u-tokyo.ac.jp/