psim言語講座（第8回）エージェントシミュレーションを書いてみる

本記事は当社が発行しているシミュレーションメールマガジンVol.8の記事です。
シミュレーションメールマガジンの詳細・購読申込はこちら のサポートページから

エージェントシミュレーション

前回は、連続シミュレーションを使って、SIR モデルを実装しました。

連続シミュレーションとは、状態変数が時間 t に対し連続的に変化するような事象を表現するものです。感染症の数理モデルの最も簡単なモデルとして、SIR モデルを紹介し、そのモデルを psim で実装してみました。

SIR モデルは、以下のような常微分方程式で表現されます。

[[ \frac{dS}{dt} = - \beta S I ]] [[ \frac{dI}{dt} = \beta S I - \gamma I ]] [[ \frac{dR}{dt} = \gamma I ]]

ここで、

• S は感受性者数
• I は感染者数
• R は回復/隔離者数

です。以下は定数です。

• $\beta$ は単位時間あたりの感染率
• $\gamma$ は単位時間あたりの回復や隔離による除去率

この式が意味しているのは、

• 感受性者は、感染者と接触すると、単位時間あたり $\beta$ の確率で感染します。接触した感染者数に比例した確率で感染します。
• 感染者は、単位時間あたり $\gamma$ の確率で回復します。

S, I, R はそれぞれの状態の人数（整数）を表していますが、SIR モデルでは、それらを連続変数としてモデリングしています。

SIR モデルのエージェントシミュレーション

今回は、より現実世界のモデルに近付け、「人」に相当するエージェントを作成し、それぞれが「状態」を持ち、適宜状態遷移していくという実装をしていきたいと思います。このようなモデリングをエージェントベースモデリングと呼びます。

まず、エージェントを定義してみます。

エージェントシミュレーションには、同期型と非同期型と呼ばれるタイプがあります。非同期型とは、この講座の最初に説明した離散イベントシミュレーションと呼ばれるタイプのシミュレーションです。非同期型はイベントを起因として状態が変化するようなシステムのシミュレーションで、イベントの発生が疎な場合に適しています。一方で同期型とは、短い固定時間ステップごとに、状態が変化するようなシステムのシミュレーションであり、イベントの発生が密な場合に適しています。

エージェントシミュレーションは、エージェント間のインタラクションが頻繁に発生する事が多いため、同期型のシミュレーションとして実装する事が多いです。SIR モデルも同期型のシミュレーションとして実装してみます。

beta = 0.001 # 単位時間あたりの感染率
gamma = 0.2 # 単位時間あたりの回復や隔離による除去率

class Agent:
    def __init__(self, state):
        # 状態は "S"/"I"/"R"
        self.state = state
    def step(self, I):
        # I は、感染者数
        if self.state == "S":
            # 感受性者は、beta * I の確率で感染する
            if np.random.uniform() <= beta * I:
                self.state = "I"
        elif self.state == "I":
            # 感染者は、gamma の確率で回復もしくは隔離
            if np.random.uniform() <= gamma:
                self.state = "R"

各エージェントは状態(self.state)を持ちます。状態は文字列 "S"/"I"/"R" の何れかとします。

同期型のエージェントシミュレーションの場合、各エージェントは必ず自身の状態を変化させる step 関数を持ちます。step 関数は短い時間間隔で呼ばれます。

step 関数の中身は、

• 感受性者は、感染者と接触すると、単位時間あたり $\beta$ の確率で感染します。接触した感染者数に比例した確率で感染します。
• 感染者は、単位時間あたり $\gamma$ の確率で回復します。

という SIR の仕様をそのまま実装しています。

次に、1000 人のエージェントを作成します。

# エージェントの初期化
agents = []

for i in range(999):
    agents.append(Agent("S"))
for i in range(1):
    agents.append(Agent("I"))
for i in range(0):
    agents.append(Agent("R"))

初期状態としては、999 人が "S" で、1 人が "I" としています。

エージェントの step 関数を定義し、初期化したら、後は時間発展させるだけです。

hist = []

def proc():
    while True:
        I = len([a for a in agents if a.state == "I"])
        # 全エージェントを状態遷移させる
        for i in range(len(agents)):
            agents[i].step(I)
        # エージェントの状態をカウントする
        S = 0
        I = 0
        R = 0
        for a in agents:
            if a.state == "S":
                S += 1
            elif a.state == "I":
                I += 1
            elif a.state == "R":
                R += 1
        # 記録
        hist.append([S, I, R])
        # 時間発展
        yield pause(1)

activate(proc)()

start(until = 30)

ここでは、proc() というプロセスを定義しています。内部に「yield pause(1)」という文がありますので、Python では proc は関数ではなく、イテレータと呼ばれます。

過去にご紹介した離散イベントシミュレーションでは、一人ごとの状態の変化をプロセスを使って表現してきました。

この例では、エージェント全員の状態変化が、ひとつのプロセスで表現されています。ステップ毎に全エージェントの状態を遷移させています。

離散イベントシミュレーションもエージェントシミュレーションも基盤として利用している仕組みは同じなのですが、実装方法が全く異なるという事がお分かり頂けたかと思います。

最後に結果をグラフ化します。

hist = np.array(hist)
plt.plot(hist[:, 0], label = "S")
plt.plot(hist[:, 1], label = "I")
plt.plot(hist[:, 2], label = "R")
plt.legend()
plt.show()

前回の連続型のシミュレーション結果と重ねてみます。

ほぼ一致する結果となっていると思います。

サンプル

エージェントシミュレーションによる SIR モデルは、こちらよりダウンロードできます。

S-Quattro Simulation System がインストールされている PC 上の適当なフォルダに上記 zip を展開下さい。

zip 内には、

• sir2.py
• sir2.bat

があります。sir2.py が SIR モデルを実装したコード例です。sir2.bat はそれを実行するための bat ファイルです。

発展

今回はエージェントシミュレーションによる SIR モデルの実装をご紹介しましたが、このモデルには様々な拡張が考えられます。

例えば、現実世界の人々は、様々な行動を取ります。例えば、ある人は朝食を家族ととり、電車に乗って職場に行き、日中はオフィスにて働き、夕刻に電車に乗って帰宅し、夜は家族と夕飯を食べます。実際には、それらの各アクションごとに、感染者数、密接状態、空調、マスク着用の割り合い、滞在時間などが異なり、感染率も変わってきます。エージェントシミュレーションでは、容易にそのような拡張が実現可能です。