ロジスティック方程式とロジスティック分布の導出　with MAXIMA

ロジスティック曲線は人口増加モデルの微分方程式から誘導され、多くの先学が各自のページで詳しく説明している。「人口増加モデル　微分方程式　ロジスティック」などと検索エンジンに打ち込めば、数学的な導出過程を説明したページが見つかる。ここでは、それらを参考にしてＭＡＸＩＭＡを使って
dy/dt=k*(1-y/L)*y 　　

という微分方程式を解きロジスティック分布を導出する。

$ \DeclareMathOperator{\abs}{abs} \newcommand{\ensuremath}[1]{\mbox{$#1$}} $

(%i2)

diffeq: 'diff(y,t)=k*y*(1-y/L);
sol1:ode2(diffeq,y,t);

\[\tag{diffeq}\frac{d}{d t} y=k y\, \left( 1-\frac{y}{L}\right) \] \[\tag{sol1}-\frac{\log{\left( y-L\right) }-\log{(y)}}{k}=t+\mathit{\% c}\]
初期値　ｔ＝０，ｙ（ｔ）＝ｃとする

(%i3)

sol2:ic1(sol1, t=0, y=c);

\[\tag{sol2}-\frac{\log{\left( y-L\right) }-\log{(y)}}{k}=\frac{k t-\log{\left( c-L\right) }+\log{(c)}}{k}\]
ｌｏｇを含む式を簡潔にすると

(%i4)

sol3:logcontract (sol2);

\[\tag{sol3}\frac{\log{\left( \frac{y}{y-L}\right) }}{k}=\frac{k t+\log{\left( \frac{c}{c-L}\right) }}{k}\]
ｙについて整理すると

(%i5)

EQ1:solve (sol3, y);

\[\tag{EQ1}[y=\frac{L c\, {{\% e}^{k t}}}{c\, {{\% e}^{k t}}-c+L}]\]

(%i6)

EQ2:rhs(EQ1[1]);;

\[\tag{EQ2}\frac{L c\, {{\% e}^{k t}}}{c\, {{\% e}^{k t}}-c+L}\]
ロジスティック分布関数を定義すると

(%i7)

define(logist(L,k,c,t),(L*c*%e^(k*t))/(c*%e^(k*t)-c+L));

\[\tag{%o7} \operatorname{logist}\left( L,k,c,t\right) :=\frac{L c\, {{\% e}^{k t}}}{c\, {{\% e}^{k t}}-c+L}\]

Ｌ＝５００，ｋ＝０．０２、初期値ｃ＝５０について
ロジスティック分布関数をグラフにすると

(%i9)

wxplot2d(logist(500,0.02,50,t),[t,-100,500]);

\[\tag{%t9} \] ロジスティック分布のグラフ

\[\tag{%o9} \]
ロジスティック分布関数をｔで微分してロジスティック確率密度関数を定義すると

(%i10)

define(logistpdf(L,k,c,t),diff(logist(L,k,c,t),t));

\[\tag{%o10} \operatorname{logistpdf}\left( L,k,c,t\right) :=\frac{L c k\, {{\% e}^{k t}}}{c\, {{\% e}^{k t}}-c+L}-\frac{L\, {{c}^{2}} k\, {{\% e}^{2 k t}}}{{{\left( c\, {{\% e}^{k t}}-c+L\right) }^{2}}}\]

Ｌ＝５００，ｋ＝０．０２、初期値ｃ＝５０について
ロジスティック確率密度関数をグラフにすると

(%i11)

wxplot2d(logistpdf(500,0.02,50,t),[t,-100,500]);

\[\tag{%t11} \] ロジスティック分布の確率密度関数のグラフ