目次
趣旨
vonBentalanffy
Gompertz
Logistic
Richards
追記
参考
趣旨
講義にて、水産資源学における色々な成長曲線を学んだので、pythonの練習も兼ねて作画してみようと思います。
なお今回はpythonの使用が目的なので、詳細な水産資源学の内容は専門書等をご覧ください。
今回使うpythonのライブラリはmatplotlib.pyplotとnumpyです。
標準で装備されているmathライブラリも使用しています。
リストの内包表現を用いて計算し、その結果を図示します。
内包表記については以下のサイトなどを併せてご覧ください。
この記事は一部、説明が丁寧でない箇所があります(プログラミングや指揮に関する詳細な記述がないです、サンプルコード程度にご高覧ください)。
気づいたら直していこうと思います()
vonBentalanffy
von Bentalanffyの成長曲線は以下の式で表される式です。
$l_{t} = l_{\infty}(1-e^{-K(t-t0)})$
$l_{t}$：t歳の体長
$l_{\infty}$：極限体長
$K$：成長係数
$t_{0}$：原点を調整するためのパラメータ
以下、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()
結果
では次に成長係数を変化させて図を描いてみます。
成長係数が0.3, 0.5, 0.7の3つの場合の曲線を比較してみます。
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()
結果
Gompertz
Gompertzの成長曲線は以下の式で表されます。
$l_{t} = l_{0}\exp(-e^{-K(t-t_{0})})$
$l_{t}$：t歳の体長
$l_{\infty}$：極限体長
$K$：成長係数
$t_{0}$：変曲点を表すパラメータ
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式に関して、y軸を対数軸に直してやるとvon Bentalanffy式の形になります。
気になる人は数式をいじってみてください。
以下ではコードと結果のみ示します。
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式に関しても複数の成長係数を比較してみましょう。
成長係数が0.3, 0.5, 0.7の3つの場合の曲線を比較してみます。
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()
結果
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で作図してみましょう、適当な変数を代入して描いてみます。
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()
結果
Gompertz式と見比べると...変曲点が少し右にずれていることに気づくかもしれません。
Logistic式の変曲点は$\frac{l_{\infty}}{2}$となることが数式から導出できます。
Logistic式についても複数の成長係数を比較してみましょう。
成長係数が0.3, 0.5, 0.7の3つの場合の曲線を比較してみます。
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()
結果
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式
- r=0の時にGompertz式
- r=1の時にLogistic式
となります。
pythonで作画してみましょう。$l_{\infty}$と$K$は適当な値を代入しておきます。
$t_{0}$は0としておきます。
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()
結果
確かにそれぞれの値で、それぞれの曲線に近しい形が出力されていることがわかります。
以上
追記
@Wolfmoon さんがnumpyを活用したコードの訂正を紹介してくださいました。
参考
専門書の例
その他