この記事では、Alex MesoudiのRで書かれたチュートリアルをpythonに翻訳する。
下記の記事ではバイアスのある・ない変異についてのモデルが紹介されている。突然変異(mutation)とは、ある形質が別の個体に継承される際に、別の形質へと変わることを意味する。

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

まずは、必要なライブラリをインストールする。

python

import pandas as pd
import numpy as np
import random
import math
import matplotlib.pyplot as plt
!pip install japanize-matplotlib
import japanize_matplotlib

バイアスのない突然変異

下記の記事で書いたコードを改変して、バイアスのない突然変異をモデル化してみる。

リスト同士の論理演算を行う場合、andを用いると後者の値をそのまま返し、orを用いると前者の値をそのまま返すというおかしな挙動をする。そのためリスト同士の論理演算を行う場合は、numpyのlogical_and()メソッドを用いる必要がある。

python

a = [True, False, False, True]
b = [True, True, False, False]

print(a and b)
#[True, True, False, False]
print(a or b)
#[True, False, False, True]
print(np.logical_and(a, b))
#[ True False False False]

python

previous_agent = agent.copy()

コードは下記のようになる。

python

def unbiased_mutation(N, mu,p_0, t_max, r_max):
    output = pd.DataFrame(np.tile(float('nan'), (t_max, r_max))) 
    output.columns = ["run_" +  str(i) for i in range(r_max)]
    for r in range(r_max):
        agent = pd.DataFrame(random.choices(["A", "B"], k=N,  weights=[p_0, 1 - p_0]), columns = ["trait"])
        output.iloc[0, r] = sum(agent["trait"] == "A") / N

        for t in range(1, t_max):
            previous_agent = agent.copy()
            mutate = np.random.rand(N)

            agent[np.logical_and(previous_agent["trait"] == "A", mutate < mu)] = "B"
            agent[np.logical_and(previous_agent["trait"] == "B", mutate < mu)] = "A"

            output.iloc[t, r] = sum(agent["trait"] == "A") / N
    
    plt.figure(figsize=(10, 5), dpi=100)
    plt.plot(output)
    plt.title(f'N = {N}')
    plt.xlabel('世代数')
    plt.ylabel("形質Aを持つエージェントの割合")
    plt.xlim([0, t_max])
    plt.ylim([0, 1])
    plt.show()

実際に、シミュレーションしてみる。p_0を0.1にすることで、形質Aを持つ個体が極めて少ない初期値にしているが、形質AとB両方が同じ割合に収束することが分かる。

python

unbiased_mutation(N=1000, mu = 0.05, p_0 = 0.1, t_max = 200, r_max = 5)

バイアスのある突然変異

続いて、バイアスのある突然変異をモデル化してみる。文化形質の場合、思い出しやすい形質へと変異することもある。
実装はAからBへ突然変異するコードを消せばよいので簡単にできる。

python

def biased_mutation(N, mu_b,p_0, t_max, r_max):
    output = pd.DataFrame([[float('nan') for i in range(t_max)] for i in range(r_max)]).T
    output.columns = ["run_" +  str(i) for i in range(r_max)]
    for r in range(r_max):
        agent = pd.DataFrame(random.choices(["A", "B"], k=N,  weights=[p_0, 1 - p_0]), columns = ["trait"])
        output.iloc[0, r] = sum(agent["trait"] == "A") / N

        for t in range(1, t_max):
            previous_agent = agent.copy()
            mutate = np.random.rand(N)

            agent[np.logical_and(previous_agent["trait"] == "B", mutate < mu_b)] = "A"

            output.iloc[t, r] = sum(agent["trait"] == "A") / N
    
    plt.plot(output)
    plt.title(f'N = {N}')
    plt.xlabel('世代数')
    plt.ylabel("形質Aを持つエージェントの割合")
    plt.xlim([0, t_max])
    plt.ylim([0, 1])
    plt.show()

実際にシミュレーションしてみよう。100世代あたりで、形質Aが集団のすべてを占めていることが見て取れる。

python

biased_mutation(N = 100, mu_b = 0.05, p_0 = 0, t_max = 200, r_max = 5)

数値シミュレーション

数値シミュレーションの場合、次のようになる。

python

p_0 = 0
t_max = 200
mu = 0.1

p = np.repeat(float('nan'), t_max)
p[0] = p_0

for i in range(1,t_max):
    p[i] = p[i-1]*(1-mu) + (1-p[i-1])*mu


plt.figure(figsize=(10, 5), dpi=100)
plt.plot(p)
plt.title(f's = {s}')
plt.xlabel('世代数')
plt.ylabel("形質Aを持つエージェントの割合")
plt.xlim([0, 150])
plt.ylim([0, 1])
plt.show()