2
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

はじめに

ロジスティック曲線は、生命科学、人口動態、化学反応、感染症モデル、機械学習の分類確率など、いろいろな場面で出てくる有名なS字曲線です。

ただ、教科書で

$$
\frac{dN}{dt} = rN\left(1-\frac{N}{K}\right)
$$

と書かれていても、最初は少し分かりにくいかもしれません。

「なぜこの式からS字になるのか?」
「指数関数的な増え方と何が違うのか?」
「$r$ や $K$ は何を意味しているのか?」
「実際のデータっぽい点に対して、どうやって曲線を当てはめるのか?」

この記事では、Pythonでダミーデータを作りながら、ロジスティック曲線を図で理解していきます。

この記事は、厳密な数理生物学の解説というより、初学者向けに「ロジスティック曲線の気持ち」をつかむための教材です。


この記事でやること

この記事では、以下を扱います。

  • ロジスティック曲線とは何か
  • 指数増殖との違い
  • 成長率 $r$ と環境収容力 $K$ の意味
  • なぜS字曲線になるのか
  • ダミーデータを作る
  • ダミーデータにロジスティック曲線をフィットする
  • 図を見ながら結果を解釈する

まず用語を整理する

成長モデルとは

成長モデルとは、ある量 $N$ が時間 $t$ とともにどのように変化するかを表すモデルです。

たとえば、$N$ は以下のような量だと思ってください。

  • 細胞数
  • 個体数
  • 菌の数
  • 売上
  • 利用者数
  • 反応生成物の量

この記事では、$N$ を「個体数」のようなものとして説明します。


指数増殖とは

一番単純な成長モデルは、指数増殖です。

指数増殖では、次のように考えます。

数が多いほど、増える量も多い。

数式では、次のように書けます。

$$
\frac{dN}{dt} = rN
$$

ここで、

  • $N$:時刻 $t$ における量
  • $t$:時間
  • $r$:成長率
  • $\frac{dN}{dt}$:$N$ の増える速さ

です。

$N$ が大きくなるほど、$rN$ も大きくなります。
そのため、指数増殖ではどんどん増えるスピードが速くなります。

しかし、現実には多くの場合、無限には増えません。

たとえば、栄養、空間、資源、競争などの制約があるからです。


ロジスティック成長とは

ロジスティック成長では、指数増殖に「増えにくさ」を加えます。

数式は次の形です。

$$
\frac{dN}{dt} = rN\left(1-\frac{N}{K}\right)
$$

ここで、新しく出てきた $K$ は 環境収容力 と呼ばれます。

環境収容力 $K$ は、ざっくり言えば、

その環境で維持できる上限のような値

です。

たとえば、$K=1000$ なら、$N$ はだいたい1000に近づいていきます。

式の中にある

$$
1-\frac{N}{K}
$$

が重要です。

この項は、$N$ が小さいときはほぼ1です。

たとえば $N \ll K$ のとき、

$$
1-\frac{N}{K} \approx 1
$$

なので、

$$
\frac{dN}{dt} \approx rN
$$

となります。

つまり、最初のうちは指数増殖に近い振る舞いをします。

一方、$N$ が $K$ に近づくと、

$$
1-\frac{N}{K} \approx 0
$$

となるので、

$$
\frac{dN}{dt} \approx 0
$$

です。

つまり、増える速さがだんだん小さくなり、最終的に頭打ちになります。


ロジスティック曲線の式

ロジスティック成長モデル

$$
\frac{dN}{dt} = rN\left(1-\frac{N}{K}\right)
$$

を解くと、次のような式になります。

$$
N(t) =
\frac{K}{1 + \left(\frac{K-N_0}{N_0}\right)e^{-rt}}
$$

ここで、

  • $N(t)$:時刻 $t$ における量
  • $N_0$:初期値、つまり $t=0$ のときの値
  • $r$:成長率
  • $K$:環境収容力
  • $e$:自然対数の底

です。

この記事では、この式を使ってロジスティック曲線を描きます。


Google Colabで実行するコード

以下のコードは、Google Colabでそのまま実行できます。

このコードでは、以下の6枚の図を作ります。

  1. 指数増殖とロジスティック成長の比較
  2. 成長速度がなぜ途中で最大になるのか
  3. 成長率 $r$ を変えたときの曲線
  4. 環境収容力 $K$ を変えたときの曲線
  5. ダミーデータへのロジスティック曲線フィット
  6. S字曲線と成長速度の関係
# ============================================
# Logistic growth tutorial
# Google Colabでそのまま実行できます
# ============================================

import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# 図の保存先
os.makedirs("fig", exist_ok=True)

# 乱数を固定して、毎回同じダミーデータが出るようにする
np.random.seed(42)

# -----------------------------
# 1. 基本設定
# -----------------------------

# 成長率
r = 0.35

# 環境収容力
K = 1000

# 初期値
N0 = 20

# 時間
t = np.linspace(0, 25, 300)


# -----------------------------
# 2. モデルの関数を定義
# -----------------------------

def exponential_growth(t, r, N0):
    """
    指数増殖モデル
    N(t) = N0 * exp(r t)
    """
    return N0 * np.exp(r * t)


def logistic_growth(t, r, K, N0):
    """
    ロジスティック成長モデル
    N(t) = K / (1 + ((K - N0)/N0) * exp(-r t))
    """
    return K / (1 + ((K - N0) / N0) * np.exp(-r * t))


# モデルの値を計算
N_exp = exponential_growth(t, r, N0)
N_log = logistic_growth(t, r, K, N0)


# -----------------------------
# 3. 図1:指数増殖とロジスティック成長
# -----------------------------

plt.figure(figsize=(8, 5))
plt.plot(t, N_exp, label="exponential growth")
plt.plot(t, N_log, label="logistic growth")
plt.axhline(K, linestyle="--", label="carrying capacity K")

plt.ylim(0, 1200)
plt.xlabel("time")
plt.ylabel("population N")
plt.title("Exponential growth vs Logistic growth")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig1_exponential_vs_logistic.png", dpi=160)
plt.show()


# -----------------------------
# 4. 図2:成長速度 dN/dt を見る
# -----------------------------

N_grid = np.linspace(0, 1200, 300)

# 指数増殖の成長速度
dNdt_exp = r * N_grid

# ロジスティック成長の成長速度
dNdt_log = r * N_grid * (1 - N_grid / K)

plt.figure(figsize=(8, 5))
plt.plot(N_grid, dNdt_exp, label="exponential: rN")
plt.plot(N_grid, dNdt_log, label="logistic: rN(1-N/K)")
plt.axvline(K, linestyle="--", label="K")
plt.axvline(K / 2, linestyle=":", label="K/2")

plt.xlabel("population N")
plt.ylabel("growth speed dN/dt")
plt.title("Why logistic growth slows down")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig2_growth_speed.png", dpi=160)
plt.show()


# -----------------------------
# 5. 図3:成長率 r を変える
# -----------------------------

plt.figure(figsize=(8, 5))

for rr in [0.15, 0.35, 0.70]:
    plt.plot(t, logistic_growth(t, rr, K, N0), label=f"r={rr}")

plt.axhline(K, linestyle="--", label="K=1000")
plt.xlabel("time")
plt.ylabel("population N")
plt.title("Effect of growth rate r")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig3_effect_of_r.png", dpi=160)
plt.show()


# -----------------------------
# 6. 図4:環境収容力 K を変える
# -----------------------------

plt.figure(figsize=(8, 5))

for KK in [500, 1000, 1500]:
    plt.plot(t, logistic_growth(t, r, KK, N0), label=f"K={KK}")

plt.xlabel("time")
plt.ylabel("population N")
plt.title("Effect of carrying capacity K")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig4_effect_of_K.png", dpi=160)
plt.show()


# -----------------------------
# 7. ダミーデータを作る
# -----------------------------

# 観測時刻
t_obs = np.linspace(0, 25, 16)

# 真のロジスティック曲線
N_true = logistic_growth(t_obs, r, K, N0)

# 観測ノイズを加える
noise = np.random.normal(0, 35, size=len(t_obs))

# 負の値にならないようにする
N_obs = np.clip(N_true + noise, 0, None)


# -----------------------------
# 8. ダミーデータにロジスティック曲線をフィットする
# -----------------------------

def logistic_for_fit(t, r, K, N0):
    """
    curve_fit用のロジスティック曲線
    """
    return K / (1 + ((K - N0) / N0) * np.exp(-r * t))


# 初期値
p0 = [0.3, 900, 30]

# パラメータの範囲
bounds = ([0.01, 100, 1], [2.0, 5000, 500])

# フィッティング
params, covariance = curve_fit(
    logistic_for_fit,
    t_obs,
    N_obs,
    p0=p0,
    bounds=bounds,
    maxfev=10000
)

r_hat, K_hat, N0_hat = params

print("===== True parameters =====")
print(f"r  = {r:.3f}")
print(f"K  = {K:.1f}")
print(f"N0 = {N0:.1f}")

print("\n===== Estimated parameters =====")
print(f"r_hat  = {r_hat:.3f}")
print(f"K_hat  = {K_hat:.1f}")
print(f"N0_hat = {N0_hat:.1f}")

# フィット曲線
N_fit = logistic_for_fit(t, r_hat, K_hat, N0_hat)


# -----------------------------
# 9. 図5:ダミーデータとフィット結果
# -----------------------------

plt.figure(figsize=(8, 5))
plt.scatter(t_obs, N_obs, label="dummy observed data")
plt.plot(t, logistic_growth(t, r, K, N0), linestyle="--", label="true curve")
plt.plot(t, N_fit, label="fitted logistic curve")

plt.xlabel("time")
plt.ylabel("population N")
plt.title("Fitting logistic curve to dummy data")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig5_dummy_data_fit.png", dpi=160)
plt.show()


# -----------------------------
# 10. 図6:S字曲線と成長速度
# -----------------------------

# ロジスティック曲線の変曲点
# N(t) = K/2 になる時刻
t_inflection = np.log((K - N0) / N0) / r

# 時間ごとの成長速度
dNdt_time = r * N_log * (1 - N_log / K)

print("\n===== Inflection point =====")
print(f"t* = {t_inflection:.2f}")
print(f"N(t*) = {logistic_growth(t_inflection, r, K, N0):.1f}")
print(f"K/2 = {K/2:.1f}")

plt.figure(figsize=(8, 5))
plt.plot(t, N_log, label="N(t)")
plt.plot(t, dNdt_time / dNdt_time.max() * K, label="scaled growth speed dN/dt")
plt.axvline(t_inflection, linestyle="--", label="inflection point")

plt.xlabel("time")
plt.ylabel("population / scaled speed")
plt.title("The middle part grows fastest")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig6_inflection_and_speed.png", dpi=160)
plt.show()

実行結果の目安

上のコードを実行すると、パラメータ推定の結果として次のような値が表示されます。

===== True parameters =====
r  = 0.350
K  = 1000.0
N0 = 20.0

===== Estimated parameters =====
r_hat  = 0.352
K_hat  = 962.3
N0_hat = 22.5

===== Inflection point =====
t* = 11.12
N(t*) = 500.0
K/2 = 500.0

図1:指数増殖とロジスティック成長の違い

fig1_exponential_vs_logistic.png

図1では、指数増殖とロジスティック成長を比較しています。

指数増殖は、時間が経つほど急激に増えていきます。
今回の設定では、指数増殖の曲線はすぐに $K=1000$ を超えてしまいます。

一方、ロジスティック成長は、最初は指数増殖に似ていますが、途中から増え方が鈍くなり、最終的に $K$ に近づいていきます。

ここで大事なのは、ロジスティック曲線は単に「上限を超えないように手で止めた曲線」ではない、ということです。

増える速さそのものが、$N$ の大きさに応じて変化します。


図2:なぜ途中で成長速度が最大になるのか

fig2_growth_speed.png

ロジスティック成長の成長速度は、

$$
\frac{dN}{dt} = rN\left(1-\frac{N}{K}\right)
$$

でした。

この式は、2つの要素に分けて考えると分かりやすいです。

まず、$rN$ は指数増殖的な要素です。

$$
rN
$$

これは、$N$ が大きいほど成長速度が大きくなる、という項です。

一方、

$$
1-\frac{N}{K}
$$

は、増えにくさを表す項です。

$N$ が小さいとき、この値は1に近いです。
そのため、最初は指数増殖に近い振る舞いになります。

しかし、$N$ が $K$ に近づくほど、この値は0に近づきます。
そのため、成長速度はだんだん小さくなります。

図2を見ると、ロジスティック成長の成長速度は、$N=K/2$ のあたりで最大になります。

今回の設定では $K=1000$ なので、

$$
\frac{K}{2} = 500
$$

のあたりで最も速く増えます。

つまり、ロジスティック曲線のS字は、

最初はゆっくり
途中で最も速く増える
最後はまたゆっくりになる

という成長速度の変化から生まれます。


少しだけ数式で見る

成長速度を

$$
g(N) = rN\left(1-\frac{N}{K}\right)
$$

とおきます。

展開すると、

$$
g(N) = rN - \frac{r}{K}N^2
$$

です。

これは $N$ に関する二次関数です。
上に凸の放物線なので、途中で最大値を持ちます。

実際に微分すると、

$$
\frac{dg}{dN} = r\left(1-\frac{2N}{K}\right)
$$

です。

成長速度が最大になるのは、

$$
\frac{dg}{dN} = 0
$$

のときなので、

$$
N = \frac{K}{2}
$$

となります。

つまり、ロジスティック成長では、量が環境収容力の半分に達したあたりで、増える速さが最大になります。


図3:成長率 $r$ を変えるとどうなるか

fig3_effect_of_r.png

図3では、$K=1000$ と $N_0=20$ は固定し、成長率 $r$ だけを変えています。

$r$ が大きいほど、曲線は早く立ち上がります。

今回の図では、

  • $r=0.15$:ゆっくり増える
  • $r=0.35$:中くらいの速さで増える
  • $r=0.70$:かなり速く増える

という違いが見えます。

ここで重要なのは、$r$ は最終的な上限を直接変えるわけではない、ということです。

この図では $K$ を固定しているので、どの曲線も最終的には $K=1000$ に近づきます。

$r$ は主に、

どれくらい速く増えるか

を決めるパラメータです。


図4:環境収容力 K を変えるとどうなるか

fig4_effect_of_K.png

図4では、成長率 $r=0.35$ と初期値 $N_0=20$ は固定し、環境収容力 $K$ だけを変えています。

$K$ が大きいほど、最終的に到達する値が大きくなります。

今回の図では、

  • $K=500$
  • $K=1000$
  • $K=1500$

を比較しています。

$K$ は、ざっくり言えば、

最終的にどのくらいの水準で頭打ちになるか

を決めるパラメータです。

ただし、実データでは $K$ の推定は簡単とは限りません。
特に、まだ頭打ちが見えていない初期のデータだけから $K$ を推定すると、不安定になりやすいです。


ダミーデータを作ってフィットしてみる

ここまでは、理想的な曲線だけを見てきました。

しかし、実際のデータにはノイズがあります。

そこで、ロジスティック曲線から人工的な観測データを作ります。

今回は、真のパラメータを

$$
r = 0.35
$$

$$
K = 1000
$$

$$
N_0 = 20
$$

としてロジスティック曲線を作り、そこにランダムなノイズを加えました。

これが、この記事で使うダミーデータです。

ダミーデータとは、実験や観測から得た本物のデータではなく、理解や検証のために人工的に作ったデータのことです。

ダミーデータの良いところは、真の答えが分かっていることです。

今回であれば、真の $r$、$K$、$N_0$ を知ったうえで、ノイズ入りデータからそれらをどれくらい推定できるかを確認できます。


図5:ダミーデータへのロジスティック曲線フィット

fig5_dummy_data_fit.png

図5では、青い点がノイズを含むダミーデータです。
破線が真のロジスティック曲線です。
実線が、ダミーデータに対してフィットしたロジスティック曲線です。

今回の実行例では、真の値と推定値は次のようになりました。

パラメータ 真の値 推定値
$r$ 0.350 0.352
$K$ 1000.0 962.3
$N_0$ 20.0 22.5

完全には一致していませんが、かなり近い値が推定されています。

特に $r$ はかなり近く推定されています。
一方で、$K$ は少し小さめに推定されています。

これは、観測ノイズがあるためです。
また、終盤のデータ点のばらつきによって、頭打ちの高さの推定は影響を受けます。

実データでは、このようなノイズや観測範囲の影響を常に考える必要があります。


図6:S字曲線と成長速度の関係

fig6_inflection_and_speed.png

図6では、ロジスティック曲線 $N(t)$ と、成長速度 $\frac{dN}{dt}$ を同じ図に描いています。

ただし、成長速度は見やすくするためにスケールを調整しています。

図を見ると、成長速度は曲線の真ん中あたりで最大になっています。

この点を 変曲点 と呼びます。

変曲点とは、曲線の曲がり方が変わる点です。

ロジスティック曲線では、変曲点は

$$
N = \frac{K}{2}
$$

のところにあります。

今回の設定では $K=1000$ なので、

$$
N = 500
$$

のところです。

コードの出力でも、

N(t*) = 500.0
K/2 = 500.0

となっています。

つまり、ロジスティック曲線では、

$N$ が $K/2$ に達したあたりで、成長速度が最大になる

という性質があります。

これがS字の中心的な理由です。


ロジスティック曲線がS字になる理由を言葉でまとめる

ロジスティック曲線がS字になる理由は、次のように説明できます。

まず、初期には $N$ が小さいです。

このとき、

$$
1-\frac{N}{K} \approx 1
$$

なので、ロジスティック成長は指数増殖に似ています。

つまり、最初は

$$
\frac{dN}{dt} \approx rN
$$

です。

そのため、$N$ が増えるほど成長速度も増えていきます。

しかし、$N$ が大きくなると、

$$
1-\frac{N}{K}
$$

が小さくなります。

この項は、資源や空間などの制約によって、増えにくくなる効果を表しています。

その結果、成長速度は途中で最大になり、その後は小さくなります。

最後には $N$ が $K$ に近づき、

$$
1-\frac{N}{K} \approx 0
$$

となるため、成長速度はほぼ0になります。

したがって、ロジスティック曲線は、

  1. 最初はゆっくり
  2. 途中で急に増える
  3. 最後は頭打ちになる

というS字になります。


指数増殖とロジスティック成長の違い

指数増殖とロジスティック成長の違いを、ざっくりまとめると次のようになります。

モデル 特徴
指数増殖 $\frac{dN}{dt}=rN$ 増えるほど、さらに速く増える
ロジスティック成長 $\frac{dN}{dt}=rN(1-\frac{N}{K})$ 最初は増えるが、やがて頭打ちになる

指数増殖は、資源や空間の制約を考えないモデルです。

一方、ロジスティック成長は、制約を非常に単純な形で取り入れたモデルです。

もちろん、現実のシステムが必ずロジスティック曲線に従うわけではありません。

しかし、

最初は増えやすいが、やがて制約によって増えにくくなる

という現象を理解するための、非常に基本的なモデルです。


実データで使うときの注意

ロジスティック曲線は便利ですが、何でもロジスティック曲線で説明できるわけではありません。

実データで使うときには、少なくとも次の点に注意が必要です。

1. 初期データだけでは K が分かりにくい

まだ頭打ちが見えていないデータだけを見ると、指数増殖なのかロジスティック成長なのか区別しにくいことがあります。

特に $K$ は、終盤のデータがないと推定が不安定になりやすいです。

2. S字に見えるからといって、同じメカニズムとは限らない

S字曲線は、いろいろな理由で現れます。

ロジスティック成長に見えるからといって、必ず

$$
\frac{dN}{dt} = rN\left(1-\frac{N}{K}\right)
$$

というメカニズムが成り立っているとは限りません。

曲線の形が似ていることと、背後の仕組みが同じであることは別です。

3. パラメータが時間とともに変わる場合もある

現実には、成長率 $r$ や環境収容力 $K$ が一定とは限りません。

環境が変われば、$r$ や $K$ も変わる可能性があります。

この記事のモデルは、$r$ と $K$ が一定であると仮定した、最も基本的なロジスティック成長モデルです。

4. ノイズや外れ値で推定結果は変わる

今回のダミーデータでも、ノイズを加えたため、推定された $K$ は真の値から少しずれました。

実データでは、測定誤差、外れ値、サンプリング間隔、観測期間などが推定結果に影響します。


おまけ:自分で試すと理解しやすい変更

上のコードを実行したあと、次の値を変えてみると理解が深まります。

成長率 r を変える

r = 0.10
r = 0.50
r = 1.00

$r$ を大きくすると、曲線の立ち上がりが急になります。


環境収容力 K を変える

K = 500
K = 2000

$K$ を大きくすると、最終的に到達する水準が大きくなります。


初期値 N_0 を変える

N0 = 5
N0 = 100
N0 = 500

$N_0$ が大きいと、最初から頭打ちに近い状態になります。


ノイズを大きくする

コード内の次の部分を探してください。

noise = np.random.normal(0, 35, size=len(t_obs))

この 35 を大きくすると、観測データのばらつきが大きくなります。

たとえば、

noise = np.random.normal(0, 100, size=len(t_obs))

にすると、フィット結果が不安定になりやすくなります。


まとめ

この記事では、ロジスティック曲線がなぜS字になるのかを、Pythonで確認しました。

ポイントは次の通りです。

  • 指数増殖は $\frac{dN}{dt}=rN$ で表される
  • ロジスティック成長は $\frac{dN}{dt}=rN(1-\frac{N}{K})$ で表される
  • $r$ は成長率、$K$ は環境収容力を表す
  • $N$ が小さいとき、ロジスティック成長は指数増殖に近い
  • $N$ が $K$ に近づくと、成長速度は0に近づく
  • 成長速度は $N=K/2$ で最大になる
  • そのため、ロジスティック曲線はS字になる
  • ダミーデータを使うと、フィッティングの挙動を確認しやすい

ロジスティック曲線は、単なる「S字の曲線」ではありません。

その背後には、

増える力と、増えにくくなる力のバランス

があります。

この視点で見ると、ロジスティック曲線の形がかなり自然に理解できると思います。

2
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?