モデルとは、事物や現象の本質的な形状や法則性を抽象化して、より単純化したものを指す言葉であり、事物や現象のモデルを作ることをモデル化という。妥当性の高いモデル化ができれば、実際に行うことが困難な実験を計算だけで行ったり、複雑な現象を再現したりするための手段としても活用できる。モデルを使って実際にどのような現象が起こるのかを予測することを一般にシミュレーションという。
モデルは、表現形式や対象の特徴によって分類することができる。静的モデルは模型や建築図面のように時間的要素を含まないもので、動的モデルはレジの待ち行列や、気象予測、生物の成長など時間的要素を含んだ事象のモデルである。また、確定モデルは不規則な現象を含まず方程式などで表せるモデルであり、確率モデルはサイコロやくじ引きのような不規則な現象を含んだモデルである。
確定モデルのひとつに複利法による預金金額の時間変化があります。
複利法は元金に利息を加算して、それを次の期間の元金として利息を計算していく方法です。
\[ (利息) = (前年の金額) \times (利率) \] \[ (今年の金額) = (前年の金額) + (利息) \]この数式モデルについて、Pythonでプログラムを書くと次のようになります。
プログラム中のリスト.append(値)は、リストの末尾に値を追加します。また、\(バックスラッシュ)は、⌥オプション＋¥(円マーク)で入力できます。
import matplotlib.pyplot as plt # matplotlib: グラフを表示するためのライブラリ
import japanize_matplotlib # japanize_matplotlib: matplotlibを日本語化するためのライブラリ
# パラメータ
rate = 0.05 # 利率
# 初期条件
amount = 10000 # 初期預金額
years = 10 # 期間
# グラフ用の配列の準備
x = [0] # リスト型変数xの0番目の要素に0
y = [amount] # リスト型変数yの0番目の要素に変数amount
# 年ごとの残高を計算
for i in range(years):
interest = int(amount * rate) # 利息 = 前年の金額 * 利率
amount = amount + interest # 今年の金額 = 前年の金額 + 利息
x.append(i + 1)
y.append(amount)
print(i + 1, "年目: \t", amount, "円") # 「\(バックスラッシュ)」は、「Option + ¥」で入力
# グラフを描画
plt.figure(figsize=(8, 5))
plt.plot(x, y, marker="o")
plt.title("複利計算のグラフ")
plt.xlabel("年数")
plt.ylabel("残高")
plt.xlim(0, years)
plt.ylim(0)
plt.grid(True)
plt.show()
1 年目: 10500 円
2 年目: 11025 円
3 年目: 11576 円
4 年目: 12154 円
5 年目: 12761 円
6 年目: 13399 円
7 年目: 14068 円
8 年目: 14771 円
9 年目: 15509 円
10 年目: 16284 円
人口16万人・１年間あたりの人口増加率5%のA市と、人口24万人・1年あたりの人口増加率1%のB市の人口を20年間シミュレーションするプログラムを作成し、何年後にB市の人口がA市の人口を上回るか求めなさい。ただし、人口増加率は変化しないものとする。
import matplotlib.pyplot as plt # matplotlib: グラフを表示するためのライブラリ
import japanize_matplotlib # japanize_matplotlib: matplotlibを日本語化するためのライブラリ
# パラメータ
rA = 0.050 # A市の人口増加率 5 %
rB = 0.010 # B市の人口増加率 1 %
# 初期条件
nA = 16 # 初期人口（万人）
nB = 24 # 初期人口（万人）
years = 20 # シミュレーション期間
# グラフ用の配列の準備
x = [0] # x軸の配列
yA = [nA] # y軸の配列（A市の人口）
yB = [nB] # y軸の配列（B市の人口）
# シミュレーションの実行
for i in range(years):
dA = rA * nA # A市の人口増加数（万人）
dB = rB * nB # B市の人口増加数（万人）
nA = nA + dA # 今年のA市の人口 = 前年のA市の人口 + A市の人口増加数 （万人）
nB = nB + dB # 今年のB市の人口 = 前年のB市の人口 + B市の人口増加数 （万人）
x.append(i + 1)
yA.append(nA)
yB.append(nB)
# グラフを描画
plt.figure(figsize=(8, 5))
plt.plot(x, yA, marker="o", label="A市")
plt.plot(x, yB, marker="o", label="B市")
plt.xlabel("時間（年）")
plt.ylabel("人口")
plt.xlim(0, years)
plt.ylim(0)
plt.grid(True)
plt.legend()
plt.show()
一次反応は、AがBに変化するときに、
を満たしながら変化する現象です。多くの自然現象や身近な現象に当てはまります。たとえば、
などが挙げられます。
単位時間あたりのAの変化率を\(k\)とすると、「Aの変化量は、そのときのAの量に比例する」ので、時間\(\Delta t\)における Aの変化量\(\Delta A\)は、 \[\Delta A = -k \times A \times \Delta t\] となります。 また、「Aの量とBの量の和が常に一定」なので、Bの変化量\(\Delta B\)は、 \[\Delta B = -\Delta A = +k \times A \times \Delta t\] となります。
import matplotlib.pyplot as plt
import japanize_matplotlib
# パラメータ
k = 0.40 # Aの変化率（1ステップあたりのAの変化量）
T = 15 # シミュレーション期間（単位は、年 / 日 / 時間 / 分 / 秒など）
dt = 0.2 # 時間間隔：1回の計算で進める時間。単位はTと同じにする。
# 初期条件
nA = 100.0 # Aの初期量
nB = 0.0 # Bの初期量
t = 0 # 時刻の初期値
# グラフ用の配列の準備
x = [t] # x軸：時刻
yA = [nA] # y軸：Aの量
yB = [nB] # y軸：Bの量
yS = [nA + nB] # y軸：Aの量 + Bの量（一定のはず）
# ループ回数
steps = int(T / dt)
# シミュレーションの実行
for i in range(steps):
dA = k * nA * dt # このステップでのAの減少量（＝Bの増加量）
nA = nA - dA # Aの量の更新
nB = nB + dA # Bの量の更新
t = t + dt # 時刻の更新
x.append(t)
yA.append(nA)
yB.append(nB)
yS.append(nA + nB)
# グラフを描画
plt.figure(figsize=(8, 5))
plt.plot(x, yA, label="A")
plt.plot(x, yB, label="B")
plt.plot(x, yS, '--', label="A + B（一定のはず）")
plt.xlabel("時間")
plt.ylabel("量")
plt.xlim(0, T)
plt.ylim(0)
plt.grid(True)
plt.legend()
plt.show()
ロジスティック成長モデルは、ある単一種の生物が一定環境内で増殖するようなときに、その生物の個体数の変動を予測するモデルです（人間の場合にも当てはまります）。何の制約もなければ、個体数は指数関数的に増加します（前項のとおり）。しかし、実際には環境や資源には限りがあるため、個体数の増加にはいずれはブレーキがかかります。つまり、個体数が増えるにつれて個体数の増加量は小さくなり、どこかの時点で飽和すると考えられます。その環境で生存しうる個体数の上限を環境収容力\(K\)と定めます。
時間\(\Delta t\)あたりの個体の変化量\(\Delta N\)は、その時点での個体数\(N\)，増加率\(r\)，環境収容力\(K\)を用いて \[\Delta N = r \times N \times (1 - \frac{N}{K}) \times \Delta t\] と表すことができます。
import matplotlib.pyplot as plt
import japanize_matplotlib
# パラメータ
r = 0.01 # 増加率（1単位時間あたりの個体数の増加割合）
K = 1000 # 環境収容力（環境が支えられる最大個体数）
T = 1000 # シミュレーション期間（単位は、年 / 日 / 時間 / 分 / 秒など）
dt = 0.5 # 時間間隔：1回の計算で進める時間。単位はTと同じにする。
# 初期条件
N = 10 # 初期個体数
t = 0 # 時刻の初期値
# グラフ用の配列の準備
x = [t] # x軸：時刻 （初期値 0）
y = [N] # y軸：個体数（初期値 N）
# ループ回数
steps = int(T / dt)
# シミュレーションの実行
for i in range(steps): # T回シミュレーションを実行
dN = r * N * (1 - N / K) * dt
N = N + dN
t = t + dt
x.append(t) # 時間のリストに新しい時間を追加
y.append(N) # 個体数のリストに新しい個体数を追加
# グラフを描画
plt.figure(figsize=(8, 5))
plt.plot(x, y)
plt.xlabel("時間")
plt.ylabel("個体数")
plt.xlim(0, T)
plt.ylim(0)
plt.title("ロジスティック成長モデル")
plt.grid(True)
plt.show()
SIRモデルは感染症の流行をシミュレーションする単純な数理モデルです。原型となったのは、W・O・カーマックとA・G・マッケンドリックによる1927年に報告された論文で、単純なSIRモデルであっても1905年のボンベイにおけるペスト流行のデータをうまく再現できることで知られています。
\(S\)，\(I\)，\(R\)の時間変化量\(\varDelta S\)，\(\varDelta I\)，\(\varDelta R\)は
\[ \varDelta S = -\beta S I \times \varDelta t \] \[ \varDelta I = (\beta S I - \gamma I) \times \varDelta t \] \[ \varDelta R = \gamma I \times \varDelta t \]と表せます。
import matplotlib.pyplot as plt
import japanize_matplotlib
# パラメータ
beta = 0.3 # 伝染率
gamma = 0.1 # 回復率
T = 100 # シミュレーション期間（日）
dt = 0.1 # 時間間隔：1回の計算で進める時間（日）。
# 初期条件
S = 0.999 # S(未感染者の割合 0 < S < 1)
I = 0.001 # I(感染者の割合 0 < I < 1)
R = 0.0 # R(回復者の割合 初期値は 0)
t = 0 # 時刻の初期値 0
# グラフ用の配列の準備
x = [0] # x軸: 時刻
yS = [S] # y軸: S(未感染者の割合)
yI = [I] # y軸: I(感染者の割合)
yR = [R] # y軸: R(回復者の割合)
# ループ回数
steps = int(T / dt)
# シミュレーションの実行
for i in range(steps):
dS = -beta * S * I * dt
dI = (beta * S * I - gamma * I) * dt
dR = gamma * I * dt
t = t + dt
S = S + dS
I = I + dI
R = R + dR
x.append(t)
yS.append(S)
yI.append(I)
yR.append(R)
# 結果のプロット
plt.figure(figsize=(8, 5))
plt.xlim(0, T)
plt.ylim(0, 1)
plt.plot(x, yS, label="$S$ 感受性保持者 (Susceptible) ")
plt.plot(x, yI, label="$I$ 感染者 (Infected)")
plt.plot(x, yR, label="$R$ 免疫保持者 (Recovered)")
plt.xlabel("時間[日]")
plt.ylabel("割合")
plt.legend()
plt.title("SIRモデル（疾病感染者数のシミュレーション）")
plt.grid(True)
plt.show()
SIRモデルは、感染症の流行の初期段階の予測、ワクチン接種の効果の評価、隔離や封じ込め策の効果のシミュレーションなど、公衆衛生の多くの場面において有用です。しかし、実際の疾病の伝播はさらに複雑であり、より高度なモデルが必要となります。