この記事では、Alex MesoudiのRで書かれたチュートリアルをpythonに翻訳する。

下記の記事ではバイアスのない伝達についてのモデルが紹介されている。これは文化進化理論におけるモデルの一つで、世代間における情報伝達の浮動(drift)をモデル化している。

モデルの概要としては、各エージェントは形質AかBのどちらかを持っており、各世代ごとにエージェントはN人いることになる。次の世代のエージェントはランダムで前の世代のエージェントの形質を継承することになるという感じだ。

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

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

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

特定の値を繰り返す行列を作る場合は、np.tile()を使う。第一引数が繰り返す値、第二引数は行、列の順に指定するタプルである。

python

pd.DataFrame(np.tile(float('nan'), (t_max, r_max)))

random.choicesでweightsによって重み付けられたリスト内の値をk回の試行回数で、ランダムに取得することができる。

python

random.choices(["A", "B"], k=N,  weights=[p_0, 1 - p_0])

下記で、unbiased_transmission関数を設計する。引数Nは個体数、p_0は形質Aの確率、t_maxは世代数、r_maxは試行回数となる。

python

#関数名が元のコードではパスカルケースとなっているが、Pythonではスネークケースにする。
def unbiased_transmission(N, 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):
        #weightsで形質を持つ確率に重み付けをする
        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()
            agent = pd.DataFrame(random.choices(previous_agent["trait"], k=N), columns = ["trait"])
            #世代における形質の割合を求め結果に格納する
            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を持つエージェントの割合")
    #matplotlibではデフォルトだと0スタートにならないため、追加で指定する必要がある。
    plt.xlim([0, t_max]) 
    plt.ylim([0, 1])
    plt.show()

では、実際にシミュレーションしてみよう。

python

unbiased_transmission(N = 10000, p_0 = 0.2, t_max = 200, r_max = 5)

このシミュレーションでの教訓は、次のように述べられている(Mesoudi, 2021)。

バイアスのない伝達それ自体は形質の頻度を変えるわけではない。
集団の規模が小さいほど、形質が消える可能性が高くなる
バイアスのない伝達は、自然・文化選択のような非ランダムなパターンを検出する上で有用である。

数値シミュレーション

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

p′=p

python

p_0 = 0.2
t_max = 200

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

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


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