psim言語講座（第9回）エージェントシミュレーションを拡張してみる

HOME
psim言語講座（第9回）エージェントシミュレーションを拡張してみる

2023年1月18日 16:36

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

SIR モデルの拡張

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

前回の記事の最後に、エージェントシミュレーションによる SIR モデルの拡張について、ご紹介しました。例えば、ある人は朝食を家族ととり、電車に乗って職場に行き、日中はオフィスにて働き、夕刻に電車に乗って帰宅し、夜は家族と夕飯を食べます。実際には、それらの各アクションごとに、感染者数、密接状態、空調、マスク着用の割り合い、滞在時間などが異なり、感染率も変わってきます。

今回は前回実装したエージェントシミュレーションに対して、そのような拡張を実装して行こうと思います。

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

まず、基盤となるエージェントのクラス定義は、前回とほとんど同じです。エージェントのグループごとに、感染率を変更できるように、beta と gamma を引数にしています。

class Agent:
    def __init__(self, state):
        # 状態は "S"/"I"/"R"
        self.state = state
    def step(self, I, beta, gamma):
        # 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"

エージェントの初期化もほぼ同じです。今回は、初期の感染者数を 10 に変更しています。

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

N = 1000

for i in range(N - 10):
    agents.append(Agent("S"))
for i in range(10):
    agents.append(Agent("I"))
for i in range(0):
    agents.append(Agent("R"))

hist = []

次にエージェントを複数のグループに分割する必要があります。そのような分割を行うためのライブラリ関数 makeGroup は以下のように実装できます。

# n 人を人数分布(dist)に従ってグループに分割する。
# 所属グループIDのリストを返す。
def makeGroup(dist, n):
    keys = np.array(list(dist.keys()))
    values = np.array(list(dist.values()))
    groups = []
    k = 0
    while len(groups) < n:
        r = np.random.uniform(0, values.sum())
        index = bisect.bisect_left(values.cumsum(), r)
        groups.extend([k] * keys[index])
        k += 1
    groups = np.array(groups[:N])
    np.random.shuffle(groups)
    return groups

makeGroup({3: 1, 4: 2}, 10)

を実行すると、例えば、

[2, 0, 2, 1, 0, 0, 0, 1, 1, 2, 2]

のような結果が得られます。これは、10人をグループ 0, 1, 2 という ID のグループに分割した結果です。リストの各要素は各自の所属グループ IDを示します。

この例では、分布を　{3: 1, 4: 2}　と指定しています、これは、10 人を 3人と 4人のグループに分割するように指示しています。その時、3人グループの数と4人グループの数の割合が 1:2 になるように分割します。ただし、確率で分割するため、正確に 1:2 になるわけではありません。

次に、グループをたばねて、各グループごとに SIR シミュレーションを行うクラスを定義し、実際にグループに分割してみます。


# 全エージェントをグループに分割したもの
# 感染はグループ内のみで発生する
class AgentGroup:
    def __init__(self, dist, beta, gamma):
        self.dist = dist
        self.beta = beta
        self.gamma = gamma
        self.ids = makeGroup(dist, N)
        self.groups = [[] for i in range(self.ids.max() + 1)]
        for i, a in enumerate(agents):
            self.groups[self.ids[i]].append(a)
    def step(self):
        S = 0
        I = 0
        R = 0
        # 各グループに対し
        for agents in self.groups:
            # グループ内の感染者数
            i = len([a for a in agents if a.state == "I"])
            # グループ内の全エージェントを状態遷移させる
            for k in range(len(agents)):
                agents[k].step(i, self.beta, self.gamma)
            # グループ内のエージェントの状態をカウントする
            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])
        return S, I, R

例えば、世帯人員分布は以下のように定義しようと思います。

世帯人数	重み
1	35.8
2	29.6
3	16.7
4	12.4
5	5.5

その場合は、以下のように定義します。

# 世帯ごとのグループ
familyGroup = AgentGroup(
    # 世帯人員分布 {世帯人員: 重み, ...}
    {1: 35.8, 2: 29.6, 3: 16.7, 4: 12.4, 5: 5.5},
    # 単位時間あたりの感染率
    0.06,
    # 単位時間あたりの回復や隔離による除去率
    0.2)

同様に、移動時のグループ、職場グループも定義します。

# 移動時のグループ
vehicleGroup = AgentGroup(
    # 移動時の乗り物の乗車人数分布 {乗車人数: 重み, ...}
    {1: 20.0, 5: 30.0, 20: 30.0, 30: 20.0},
    # 単位時間あたりの感染率
    0.04,
    # 単位時間あたりの回復や隔離による除去率
    0.2)
# 職場グループ
officeGroup = AgentGroup(
    # 日中の職場(学校)の人数 {乗車人数: 重み, ...}
    {1: 20.0, 5: 30.0, 20: 30.0, 30: 20.0},
    # 単位時間あたりの感染率
    0.03,
    # 単位時間あたりの回復や隔離による除去率
    0.2)

これらを統合して最終的な SIR シミュレーションでは、

世帯ごとのグループで 4 回
移動時のグループで 1 回
職場グループを 4 回
移動時のグループで 1 回

の状態遷移を繰り返します。

def proc():
    while True:
        # 世帯での感染
        for k in range(4):
            S, I, R = familyGroup.step()
            yield pause(1)
        # 出社移動時の感染
        for k in range(1):
            S, I, R = vehicleGroup.step()
            yield pause(1)
        # 日中の職場での感染
        for k in range(4):
            S, I, R = officeGroup.step()
            yield pause(1)
        # 帰宅移動時の感染
        for k in range(1):
            S, I, R = vehicleGroup.step()
            yield pause(1)

activate(proc)()

start(until = 100)

実行結果は以下のようになりました。