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?

Day 6

NSSOL Advent Calendar 2025

@yamamorisobain日鉄ソリューションズ株式会社

【続編】ロトカ・ヴォルテラを確率過程で再現 ― Gillespie法＋ランダムテレポート捕食モデル

2
Last updated at Posted at 2025-12-05

【続編】ロトカ・ヴォルテラ方程式 を “確率過程 & 個体ベース” で再現してみる

― Gillespie 法・自律MOB・ランダムテレポート捕食モデルまで ―

0. はじめに

この記事は、NSSOL Advent Calendar6日目の記事で、2日目の記事の続編です。
前回の記事では、Minecraft Education Edition ＋ MakeCode を使って、
ロトカ・ヴォルテラ 方程式（捕食者–被食者モデル）を数値的に解き、可視化しました。

しかし、前回モデルには次のような「限界」があります：

  • 個体数が実数（1.2体など）
  • イベント（捕食・出生）が連続時間
  • 時間発展は決定論的
  • MOBの増減は、数値計算の結果を反映する形

せっかくMinecraftで可視化するのであれば、「もっと“Minecraft の世界”っぽい、個体（MOB）としてのランダムな振舞いや、MOB間の相互作用を反映した系にしたい」という欲が出てくるわけです。

そこで本記事では：

  • ロトカ・ヴォルテラ方程式で描かれる捕食・被食のふるまいを 確率過程（Gillespie法） にし、
  • MOB の数が整数として揺らぐモデル を考え、その統計的なふるまいを調べました
  • さらに発展として 自律行動する捕食者（ランダムテレポート＋局所捕食） を検討します

※自律行動する捕食者部分は検討までで、実装は間に合いませんでした…

結論としては、Minecraft（教育版） でも ロトカ・ヴォルテラ にかなり近い「確率ゆらぎのある生態系」を作れそうだ、ということが分かりました。

1. 前回モデル（連続 LV）の復習

ロトカ・ヴォルテラ 方程式（標準形）

$$
\begin{aligned}
\frac{dx}{dt} &= a x - b xy \\
\frac{dy}{dt} &= c xy - d y
\end{aligned}
$$

  • 被食者 x は、時間発展で増え（ax項）、捕食によって減少する（bxy項）
  • 捕食者 y は、捕食によって増え(cxy項)、生存コストで減る（dy項）
  • 時間とともに「振動（サイクル）」が現れる

前回はこの式をオイラー法で解いた結果を離散化し、マインクラフトのMobの増減として表現する方法を紹介しました。

その結果、被食者と捕食者の数が周期的に変動する様子を確認できましたが、一方で、

  • 空間構造なし
  • 個体ゆらぎなし
  • 常に滑らかな軌道

という「数理モデルとしては正しいが、生態系としては物足りない」側面があります。

2. ロトカ・ヴォルテラ を確率過程化 ― Gillespie 法とは？

実世界では、出生・捕食・死亡は連続ではなく離散イベントです。連続時間の方程式を、離散イベントの確率過程として扱うため、化学反応や生態系モデルで広く使われるGillespie法を採用しました

基本イベントは3つ

Gillespie法でのシミュレーションにおいて、今回想定するイベントは以下の3つです。
※捕食による被食者の減少と捕食者の増加は、ロトカ・ヴォルテラ方程式では、-bxyとcxyのように、異なる定数で規定されていますが、今回はK3Tunnelのサンプル内で採用されているb=c=2を前提にシミュレーションを行うため、同一のイベントとして処理しています。

イベント 反応式 レート 意味
① 被食者の出生 X → 2X a·X 被食者が生まれる
② 捕食 X+Y → 2Y b·X·Y 捕食により X が減り Y が増える※
③ 捕食者の死 Y → ∅ d·Y 捕食者が死ぬ

Gillespie 法のステップ

Gillespie法でのシミュレーションは以下のステップで行います。

  1. その時点のレート（発生しやすさ） λ₁, λ₂, λ₃ を計算
  2. すべてのイベントの発生確率の合計 λ = λ₁+λ₂+λ₃
  3. 乱数を使って「次に起こるイベント」を1つ選ぶ
  4. イベントを実行（X, Y を更新）
  5. 1～4を繰り返す

このステップの特徴は、
→ 1度に起こるのは1イベントだけ
→ 発生するイベントが確率的に選ばれるため、個体数のゆらぎや絶滅が発生しうる
という点です。このアプローチで、ある初期条件からの個体群の時間発展を複数回シミュレーションし、各時間帯の個体数を統計的に処理することで、その系の時間発展の振舞いをうかがい知ることができるようになります。

3. python で Gillespie 法を実装する

Gillespie法でのシミュレーションをpythonで実装します。ロトカ・ヴォルテラ方程式の数値解と比較するできるように、ロトカ・ヴォルテラ方程式も併せて解きます（前回はオイラー法でしたが、今回はルンゲ・クッタ法を使っています）
なお、このコードは、ChatGPTで生成したものをClaudeでレビューしています。

シミュレーション用コードサンプル（長いので折りたたみ）

コードを開く 
import random
import math
import csv


# ============================
# 1. モデルパラメータ
# ============================
# ロトカ・ヴォルテラ方程式：
#   dx/dt = a x - b x y
#   dy/dt = c x y - d y
#
# ここで x, y は「密度」。
# 実際の個体数は N = K x, P = K y とする。

# ★ユーザー指定の係数
A = 1.0   # a
B = 2.0   # b
C = 2.0   # c
D = 0.8   # d

# 系のサイズ（個体数スケーリング）
# 必要に応じて K を大きく/小さくして調整してください
K = 10000  # 例: 1000。大きくすると揺らぎが小さくなるが時間もかかる。

# ★ユーザー指定の初期密度
x0 = 1.0   # 初期獲物密度
y0 = 0.2   # 初期捕食者密度

# 整数個体数に変換
N0 = int(K * x0)
P0 = int(K * y0)

# シミュレーション条件
T_MAX = 20.0      # 観測する時間の上限
DT_SAMPLE = 0.1   # この刻みで「密度」を記録

# アンサンブル（試行回数）
N_TRAJ = 100      # 多いほど平均が滑らかになるが計算時間↑


# ============================
# 2. Gillespie 1 本分のシミュレーション
# ============================
def simulate_one_gillespie(seed=None):
    """
    スケーリング付き Gillespie 法による 1 本分の軌道。
    密度 x=N/K, y=P/K を、DT_SAMPLE ごとの等間隔時刻で返す。
    """
    if seed is not None:
        random.seed(seed)

    t = 0.0
    N = N0
    P = P0

    # サンプル時刻リスト
    times = []
    t_sample = 0.0
    while t_sample <= T_MAX + 1e-9:
        times.append(t_sample)
        t_sample += DT_SAMPLE

    xs = []  # 密度 x = N/K
    ys = []  # 密度 y = P/K

    sample_index = 0

    while t < T_MAX and (N > 0 or P > 0):
        # 反応率（スケーリング付き）
        # N = K x, P = K y とすると、
        #   捕食関連の反応率は (b/K) N P, (c/K) N P になる。
        r1 = A * N             # 獲物出生
        r2 = (B / K) * N * P   # 捕食による獲物減少
        r3 = (C / K) * N * P   # 捕食者出生
        r4 = D * P             # 捕食者死亡

        R = r1 + r2 + r3 + r4
        if R <= 0.0:
            break  # これ以上何も起きない

        u1 = 1.0 - random.random()
        tau = -math.log(u1) / R

        # tau の間に通過するサンプル時刻を、現在の状態で埋める
        while sample_index < len(times) and times[sample_index] <= t + tau:
            xs.append(N / K)
            ys.append(P / K)
            sample_index += 1

        t += tau
        if t > T_MAX:
            break

        # どのイベントが起きるか
        u2 = random.random() * R
        if u2 < r1:
            # 獲物出生
            N += 1
        elif u2 < r1 + r2:
            # 捕食 → 獲物減少
            if N > 0:
                N -= 1
        elif u2 < r1 + r2 + r3:
            # 捕食者出生
            P += 1
        else:
            # 捕食者死亡
            if P > 0:
                P -= 1

    # シミュレーションが終わった後、残りのサンプル時刻を最後の状態で埋める
    while sample_index < len(times):
        xs.append(N / K)
        ys.append(P / K)
        sample_index += 1

    return times, xs, ys


# ============================
# 3. アンサンブル平均・標準偏差 + 全軌道保存
# ============================
def simulate_ensemble():
    """
    Gillespie を N_TRAJ 本走らせて、
      - base_times: サンプル時刻
      - mean_x, mean_y: 各時刻の密度の平均
      - std_x, std_y: 各時刻の密度の標準偏差
      - trajectories_x, trajectories_y:
            各試行の x(t), y(t) を格納した2次元リスト
    を返す。
    """
    print(f"Gillespie アンサンブル計算開始: 本数 = {N_TRAJ}")

    # まず 1本走らせて時間軸を確定
    base_times, xs0, ys0 = simulate_one_gillespie(seed=0)
    n_samples = len(base_times)

    # 全軌道を保存するリスト
    trajectories_x = [xs0]
    trajectories_y = [ys0]

    # 平均・分散用
    sum_x = [0.0] * n_samples
    sum_y = [0.0] * n_samples
    sum_x2 = [0.0] * n_samples
    sum_y2 = [0.0] * n_samples

    # 1本目を反映
    for i in range(n_samples):
        x = xs0[i]
        y = ys0[i]
        sum_x[i] += x
        sum_y[i] += y
        sum_x2[i] += x * x
        sum_y2[i] += y * y

    # 残り N_TRAJ-1 本を計算
    for k in range(1, N_TRAJ):
        times_k, xs_k, ys_k = simulate_one_gillespie(seed=k)
        trajectories_x.append(xs_k)
        trajectories_y.append(ys_k)

        for i in range(n_samples):
            x = xs_k[i]
            y = ys_k[i]
            sum_x[i] += x
            sum_y[i] += y
            sum_x2[i] += x * x
            sum_y2[i] += y * y

        # ★進捗表示：何本目まで終わったか
        done = k + 1
        percent = 100.0 * done / N_TRAJ
        # 必要に応じて頻度を絞りたい場合は if done % 5 == 0: などで絞ってもよい
        print(f"  進捗: {done}/{N_TRAJ} 本 ({percent:.1f}%) 完了")

    # 平均と標準偏差を計算
    mean_x = []
    mean_y = []
    std_x = []
    std_y = []

    for i in range(n_samples):
        mx = sum_x[i] / N_TRAJ
        my = sum_y[i] / N_TRAJ
        vx = sum_x2[i] / N_TRAJ - mx * mx
        vy = sum_y2[i] / N_TRAJ - my * my
        if vx < 0:
            vx = 0.0
        if vy < 0:
            vy = 0.0
        mean_x.append(mx)
        mean_y.append(my)
        std_x.append(math.sqrt(vx))
        std_y.append(math.sqrt(vy))

    print("Gillespie アンサンブル計算完了。")
    return base_times, mean_x, mean_y, std_x, std_y, trajectories_x, trajectories_y


# ============================
# 4. ロトカ・ヴォルテラ方程式（ODE）の数値解
# ============================
def solve_lv_ode():
    """
    ロトカ・ヴォルテラ方程式を密度 x,y について解く。
    サンプル時刻と同じ刻みで x,y を返す。
    （単純な RK4 で DT_SAMPLE 刻み）
    """
    print("LV ODE を解いています...")

    times = []
    xs = []
    ys = []

    t = 0.0
    x = x0
    y = y0

    while t <= T_MAX + 1e-9:
        times.append(t)
        xs.append(x)
        ys.append(y)

        # 4次のRunge-Kutta法（1ステップ DT_SAMPLE）
        h = DT_SAMPLE

        def f_x(x_val, y_val):
            return A * x_val - B * x_val * y_val

        def f_y(x_val, y_val):
            return C * x_val * y_val - D * y_val

        k1x = f_x(x, y)
        k1y = f_y(x, y)

        k2x = f_x(x + 0.5 * h * k1x, y + 0.5 * h * k1y)
        k2y = f_y(x + 0.5 * h * k1x, y + 0.5 * h * k1y)

        k3x = f_x(x + 0.5 * h * k2x, y + 0.5 * h * k2y)
        k3y = f_y(x + 0.5 * h * k2x, y + 0.5 * h * k2y)

        k4x = f_x(x + h * k3x, y + h * k3y)
        k4y = f_y(x + h * k3x, y + h * k3y)

        x = x + (h / 6.0) * (k1x + 2 * k2x + 2 * k3x + k4x)
        y = y + (h / 6.0) * (k1y + 2 * k2y + 2 * k3y + k4y)
        t += h

    print("LV ODE 計算完了。")
    return times, xs, ys


# ============================
# 5. メイン：実行・CSV 出力・プロット
# ============================
def main():
    # Gillespie アンサンブル
    times, mean_x, mean_y, std_x, std_y, trajectories_x, trajectories_y = simulate_ensemble()

    # ODE
    ode_times, ode_x, ode_y = solve_lv_ode()

    # --- アンサンブル結果 + ODE を CSV 出力 ---
    filename_ensemble = "lv_gillespie_scaled.csv"
    with open(filename_ensemble, "w", newline="") as f:
        writer = csv.writer(f)
        writer.writerow([
            "time",
            "ode_x", "ode_y",
            "mean_x", "mean_y",
            "std_x", "std_y"
        ])
        for t, ox, oy, mx, my, sx, sy in zip(times, ode_x, ode_y, mean_x, mean_y, std_x, std_y):
            writer.writerow([t, ox, oy, mx, my, sx, sy])

    print("アンサンブル平均 + ODE の CSV 出力完了:", filename_ensemble)

    # --- 各試行の軌道を別 CSV に出力 ---
    filename_traj = "lv_gillespie_trajectories.csv"
    with open(filename_traj, "w", newline="") as f:
        writer = csv.writer(f)
        writer.writerow(["traj_index", "time", "x", "y"])
        for k in range(N_TRAJ):
            xs_k = trajectories_x[k]
            ys_k = trajectories_y[k]
            for i, t in enumerate(times):
                writer.writerow([k, t, xs_k[i], ys_k[i]])

    print("各試行の軌道 CSV 出力完了:", filename_traj)

    print("  パラメータ: a=%.3f, b=%.3f, c=%.3f, d=%.3f, K=%d" % (A, B, C, D, K))
    print("  初期密度: x0=%.3f, y0=%.3f  => N0=%d, P0=%d" % (x0, y0, N0, P0))
    print("  試行回数: N_TRAJ=%d" % N_TRAJ)
    print("  時間: 0 ～ %.2f, サンプル刻み=%.3f" % (T_MAX, DT_SAMPLE))

if __name__ == "__main__":
    main()

4. 実行結果：ゆらぎ・絶滅・“個体らしさ” が見えてくる

このプログラムの実行結果について説明します。まず、乱数を用いて確率的に時間発展するモデルのため、特に個体数が少ないケースでは、揺らぎの影響が大きくロトカ・ヴォルテラ方程式のふるまいから大きく外れることがあります。極端なケースとして、捕食者や被食者が絶滅した試行の例をプロットします。（以降の示すグラフでは、横軸が時間、縦軸が個体密度です）

①オオカミが絶滅して、ヒツジが爆増
t=17あたりでオオカミが絶滅しています。捕食者がいなくなった結果、この後ヒツジは爆増します。
image.png

②ヒツジが絶滅して、オオカミも絶滅
t=13あたりでヒツジが食い尽くされて絶滅しました。結果、後を追うようにオオカミも絶滅しています。
image.png

今回の取り組みでは、同じ個体数での試行に関して100回のシミュレーション結果の平均をとることで揺らぎの影響を抑えています。このモデルは、個体数が無限大となる極限でロトカ・ヴォルテラ方程式の解に一致することが期待値なので、シミュレーションする群れの個体数を大きくしていき、その結果をロトカ・ヴォルテラ方程式の数値計算結果と比較します。

N=100 P=20

シミュレーションの後半で、ヒツジの数が発散するケースの影響が大きく出ている。

被食者 捕食者
x100.png y100.png

N=1000 P=200

1周期目は、ロトカ・ヴォルテラ方程式の解にほぼ一致、2周期目以降で乖離が大きくなっている

被食者 捕食者
x1000.png y1000.png

N=10000 P=2000

後半の周期でずれが見られるが、個体数が小さいケースよりも乖離が小さい

被食者 捕食者
x10000.png y10000.png

N=100000 P=20000

シミュレーションした時間の範囲では、ロトカ・ヴォルテラ方程式の解にほぼ一致

被食者 捕食者
x100000.png y100000.png

5. Minecraftの世界でロトカ・ヴォルテラに近い振舞いを再現するには？

→ ランダムテレポート × 局所捕食 というアイデア

これまで、扱ってきたロトカ・ヴォルテラ方程式や、Gillespieモデルでのシミュレーションでは、捕食者と被食者の相互作用項がbxyのように単純な積となっています。これは捕食者と被食者が互いの位置関係によらず一様に相互作用する完全混合の状態を前提としていることを意味しています。

一方で Minecraft の世界は、

  • 捕食者は「近くの被食者」しか捕食（攻撃）しない
  • 捕食者と被食者が遠く離れていると遭遇しない
  • 地形が影響する

というように、 完全混合（ロトカ・ヴォルテラの前提）ではない 世界観となっています。ロトカ・ヴォルテラ方程式の振舞いに近い系を作ろうとすると、すべてのMOB同士が距離などの位置関係に依存せず相互作用する状況を整える必要があります。

そこで有望と考えられるアプローチが「捕食者をワールド内でランダムテレポートさせる」手法です。

  1. 捕食者（例：wolf）をランダムな位置にテレポート
  2. 半径R以内に被食者（sheep）がいるか調べる
  3. いれば確率 p で kill → 捕食成功
  4. 捕食成功したら一定確率で増殖

この方法は、以下の点でロトカ・ヴォルテラの完全混合を近似できる可能性があります。

  • 捕食者が「どこでも均等に羊を探索する」
  • → 遭遇率が x×y に近づく（ロトカ・ヴォルテラ方程式の前提に近づく）
  • （テレポートの頻度を適切に高めることで）高速でランダムに混ざるため、空間の偏りが弱くなる

6. 自律MOBモデル（簡易版）MakeCodeコード例（未完）

オオカミがランダムテレポートし、近傍のヒツジを一定確率で捕食するパターンをサンプルコードとして提示したかったのですが、今回はタイムアップでそこまで至りませんでした。MakeCodeでの実装例と実行結果については、続々編として、後日記事にまとめたいと思います。

7. まとめ

  • Gillespie 法で「確率的 ロトカ・ヴォルテラ」を実装した
  • 生態系らしいゆらぎ・絶滅が自然に表れた
  • Minecraft でも “完全混合” をランダムテレポートで近似できる可能性を検討した
  • MakeCodeでの実装については未完
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
Sign upLogin
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?