微分方程式のパラメータフィッティング

出会い

ある日、面白い本に出会った。

微分方程式で数学モデルを作ろう

ある日、たまたま行ったコワーキングスペース 5.6に本棚があった。そこにある本を何気なく取ってパラパラとめくったところ、ビビビッとくるものがあり、即購入した一冊である。

特に感銘を受けたのが以下の図と文章である。

図の一番左の列は現実の世界を、一番右の列は数学の世界をあらわし、これら２つの世界を結びつけるものが中央の列にある。応用数学では、時に１，２，３の枠をお座なりに通過することもあるが、ほとんどの数学教育は４番目の枠にだけかかわっていて、いずれにしても５，６，７の枠へと進むことはない。数学教育の中で、この１から７までのつながりを重要視しない限り、現代の科学・工業社会のなかで当面している多くの問題の解決に数学が重要な役割を果しえていることを、学生は十分に理解できないだろう。

そう！！本当にそう！！「ぶっちゃけ数学使わなくね？」と思うのは、実世界を数学にモデリングする訓練をしないし、それを改善するという営みをしないからだ。数学教育では数式の操作に終始し、使い方のよくわからない道具を勉強させられている感が自分には強かった。この本では、いろいろな現象を微分方程式の方面から観察し、数式を改善したり、数学的な操作で解を求めていく一冊となっていて、個人的には膝を打つ一冊であった。

ヴェアフルストの人口モデル

書籍には、一例としてヴェアフルストの人口モデルというものが記載されていた。
ヴェアフルストの人口モデルは、人口の増加が環境資源によって制約される現象を数学的に記述したものである。

\frac{dN}{dt} = rN \left( 1 - \frac{N}{C} \right)

\begin{align}
&N: \text{時間 } t \text{ における人口サイズ} \\
&r: \text{成長率} \\
&C: \text{環境収容力（人口の上限を決定するパラメータ）} \\
&\frac{dN}{dt}: \text{人口の時間に対する変化率}\ \\
\end{align}

そのあと、書籍ではデータが示され、

年度	人口 (×10^6)
1820	9.6
1830	12.9
1840	17.1
1850	23.2
1860	31.4
1870	38.6
1880	50.2
1890	62.9
1900	76.0
1910	92.0
1920	106.5
1930	123.2

フィッティングしたグラフとパラメーターが示されている。

このヴェアフルストの人口モデルは解析解が分かっている。そのため、その解析解に対してフィッティングは可能である。

N = \frac{N_{\infty}}{1 + \left[ \left( \frac{N_{\infty}}{N_0} \right) - 1 \right] e^{-rt}}

では、**解析解のわからない微分方程式のパラメータフィッティングってどうやってやるんだ？**というところが、今日の記事の内容である。
今回の記事であるが、おおむねChatGPTに解法のお伺いを立て、それをまとめた内容となっている。

疑似中央差分によるフィッティング

ヴェアフルストの人口モデルの式を考えるとき、問題となるのは微分である。ここをどう扱うかという問題に帰着する。

\frac{dN}{dt} = rN \left( 1 - \frac{N}{C} \right)

最初に私が思いついたのは中央差分である。

ある注目するi番目の値の微分は、i-1番目とi+1番目を結ぶ直線と同じ傾きをしていると仮定する。

\frac{dN[i]}{dt} \simeq \frac{N[i+1]-N[i-1]}{t[i+1]-t[i-1]}

これを用いて、ChatGPTさんにコードを出してもらった。私はほとんど手を動かしておらず、プロンプト芸であーだこーだ言ってるだけである。

疑似中央差分によるフィッティングのコード

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# 必要なライブラリのインポート
!apt-get -y install fonts-ipafont-gothic
import matplotlib

# フォント設定（日本語は不要）
plt.rcParams['axes.unicode_minus'] = False

# 年度と人口データ
years = np.array([1820, 1830, 1840, 1850, 1860, 1870, 1880, 1890, 1900, 1910, 1920, 1930])
population = np.array([9.6, 12.9, 17.1, 23.2, 31.4, 38.6, 50.2, 62.9, 76.0, 92.0, 106.5, 123.2])

# データを前半と後半に分割
split_index = len(years) // 2
years_train, years_test = years[:split_index], years[split_index:]
population_train, population_test = population[:split_index], population[split_index:]

# 差分近似で微分値を計算する関数
def calculate_derivatives(years, population):
    dN_dt = []
    for i in range(1, len(population) - 1):
        derivative = (population[i + 1] - population[i - 1]) / (years[i + 1] - years[i - 1])
        dN_dt.append(derivative)
    return np.array(dN_dt)

# 微分値と中間のNデータ
dN_dt_train = calculate_derivatives(years_train, population_train)
N_mid_train = population_train[1:-1]  # 中間の人口値

# 損失関数の定義
def loss_function(params):
    r, C = params
    model_derivative = r * N_mid_train * (1 - N_mid_train / C)
    return np.sum((dN_dt_train - model_derivative) ** 2)

# 初期値と境界条件
initial_guess = [0.02, 200]  # rとCの初期値
bounds = [(0, None), (0, None)]  # rとCは正の値のみ

# 最適化
result = minimize(loss_function, initial_guess, bounds=bounds)

# 最適化結果
r_opt, C_opt = result.x
min_loss = result.fun  # 最小化された損失関数の値
print(f"Estimated r: {r_opt}")
print(f"Estimated C: {C_opt}")
print(f"Minimized loss (error metric): {min_loss}")

# 推定されたモデルを差分方程式でシミュレーション
def simulate_population_difference(t, N0, r, C):
    N_sim = [N0]
    dt = t[1] - t[0]
    for i in range(len(t) - 1):
        N_current = N_sim[-1]
        dN = r * N_current * (1 - N_current / C) * dt
        N_sim.append(N_current + dN)
    return np.array(N_sim)

# 全データに対するシミュレーション
predicted_population = simulate_population_difference(years, population_train[0], r_opt, C_opt)

# グラフの表示
plt.figure(figsize=(10, 6))
plt.plot(years, population, 'o', label="Observed Data", color='blue')
plt.plot(years, predicted_population, '-', label="Simulation (Entire Data)", color='orange')

# データの分割点に縦線を描画
plt.axvline(x=years[split_index], color='red', linestyle='--', label="Train-Test Split")

# ラベルとタイトルの設定
plt.xlabel("Year")
plt.ylabel("Population (×10^6)")
plt.title("Population Prediction: Observed vs Simulation")
plt.legend()
plt.grid()
plt.show()

Estimated r: 0.034560641446984744
Estimated C: 129.2386726853784
Minimized loss (error metric): 0.006618323987894751

全体を前後ろに分け、前半部分でパラメータのフィッティングを行い、差分法でシミュレーションを行っている。どうやらこれでもぱらっと出るらしい。ChatGPTすごい。

平均微分によるフィッティング

次に思いついたのが微分の平均である。

注目しているi番目の値に対して、i-1番目とi番目の傾きを計算する。そのあと、i番目とi+1番目の傾きを計算する。それらの平均値を取る。

\frac{dN[i]}{dt} \simeq \frac{1}{2} \lbrace \frac{N[i]-N[i-1]}{t[i]-t[i-1]} + \frac{N[i+1]-N[i]}{t[i+1]-t[i]} \rbrace

この値を使って計算してみる。

平均微分によるフィッティングのコード

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# 必要なライブラリのインポート
!apt-get -y install fonts-ipafont-gothic
import matplotlib

# フォント設定（日本語は不要）
plt.rcParams['axes.unicode_minus'] = False

# 年度と人口データ
years = np.array([1820, 1830, 1840, 1850, 1860, 1870, 1880, 1890, 1900, 1910, 1920, 1930])
population = np.array([9.6, 12.9, 17.1, 23.2, 31.4, 38.6, 50.2, 62.9, 76.0, 92.0, 106.5, 123.2])

# データを前半と後半に分割
split_index = len(years) // 2
years_train, years_test = years[:split_index], years[split_index:]
population_train, population_test = population[:split_index], population[split_index:]

# 差分近似で微分値を計算する新しい関数
def calculate_derivatives_new(years, population):
    dN_dt = []
    for i in range(1, len(population) - 1):
        # 新しい近似式を適用
        derivative = 0.5 * (
            (population[i + 1] - population[i]) / (years[i + 1] - years[i]) + 
            (population[i] - population[i - 1]) / (years[i] - years[i - 1])
        )
        dN_dt.append(derivative)
    return np.array(dN_dt)

# 新しい微分値と中間のNデータ
dN_dt_train_new = calculate_derivatives_new(years_train, population_train)
N_mid_train_new = population_train[1:-1]  # 中間の人口値

# 損失関数の定義（新しい近似式を使用）
def loss_function_new(params):
    r, C = params
    model_derivative = r * N_mid_train_new * (1 - N_mid_train_new / C)
    return np.sum((dN_dt_train_new - model_derivative) ** 2)

# 初期値と境界条件
initial_guess = [0.02, 200]  # rとCの初期値
bounds = [(0, None), (0, None)]  # rとCは正の値のみ

# 最適化
result_new = minimize(loss_function_new, initial_guess, bounds=bounds)

# 最適化結果
r_opt_new, C_opt_new = result_new.x
min_loss_new = result_new.fun  # 最小化された損失関数の値
print(f"Estimated r (new): {r_opt_new}")
print(f"Estimated C (new): {C_opt_new}")
print(f"Minimized loss (new): {min_loss_new}")

# 推定されたモデルを差分方程式でシミュレーション
def simulate_population_difference(t, N0, r, C):
    N_sim = [N0]
    dt = t[1] - t[0]
    for i in range(len(t) - 1):
        N_current = N_sim[-1]
        dN = r * N_current * (1 - N_current / C) * dt
        N_sim.append(N_current + dN)
    return np.array(N_sim)

# 前半のデータを使用してシミュレーション
predicted_population = simulate_population_difference(years, population_train[0], r_opt_new, C_opt_new)

# 逐次更新によるシミュレーション
def sequential_simulation(t, N0, r, C, observed_data):
    N_sim = [N0]
    for i in range(1, len(t)):
        dt = t[i] - t[i - 1]
        N_current = observed_data[i - 1]
        # 観測データを使用してdN/dtを更新
        dN_dt = r * N_current * (1 - N_current / C)
        N_next = N_current + dN_dt * dt
        N_sim.append(N_next)
    return np.array(N_sim)

# 逐次更新シミュレーションを実行
sequential_population = sequential_simulation(years, population_train[0], r_opt_new, C_opt_new, population)

# グラフの表示
plt.figure(figsize=(10, 6))

# 観測データとシミュレーションの結果
plt.plot(years, population, 'o', label="Observed Data", color='blue')
plt.plot(years, predicted_population, '-', label="Simulation (Entire Data)", color='orange')
#plt.plot(years, sequential_population, '--', label="Sequential Simulation", color='green')

# データの分割点に縦線を描画
plt.axvline(x=years[split_index], color='red', linestyle='--', label="Train-Test Split")

# ラベルとタイトルの設定
plt.xlabel("Year")
plt.ylabel("Population (×10^6)")
plt.title("Population Prediction: Observed vs Simulation")
plt.legend()
plt.grid()
plt.show()

# 誤差を評価する（後半のデータに対する誤差）
error = np.abs(population_test - sequential_population[split_index:])
error_metric = np.mean(error)  # 平均絶対誤差

print(f"Error Metric (Mean Absolute Error) for Sequential Simulation: {error_metric:.2f}")

Estimated r (new): 0.034560641446984744
Estimated C (new): 129.2386726853784
Minimized loss (new): 0.006618323987894751

計算してみた結果、値は変わらなかった。

ノイズに対する耐性

ChatGPTに色々と相談をしながらプログラムを組んできた。すると気になることを言っている。

ノイズがあるとこの方法はだめらしい。では、データを水増ししてノイズを入れて計算してみよう。
プロンプトはかなり雑でもなんとかしてくれた。

年度と人口データを水増ししたいです。3次スプライン補完を行うことで、1年刻みのデータに拡張してください。

スプライン補完したpopulationのデータにN(0,5)の白色ノイズを加えてシミュレーションするコードを書いて下さい。

というわけで三次スプライン補完をした後、ノイズを入れたデータにフィッティングをしてみる。

ノイズ耐性の検証コード

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.interpolate import CubicSpline

# 年度と人口データ（元データ）
years = np.array([1820, 1830, 1840, 1850, 1860, 1870, 1880, 1890, 1900, 1910, 1920, 1930])
population = np.array([9.6, 12.9, 17.1, 23.2, 31.4, 38.6, 50.2, 62.9, 76.0, 92.0, 106.5, 123.2])

# 3次スプライン補完を行う
spline = CubicSpline(years, population)

# 新しい1年刻みの年度データを作成
extended_years = np.arange(years[0], years[-1] + 1, 1)

# 1年刻みで補完した人口データを計算
extended_population = spline(extended_years)

# N(0,5)の白色ノイズを加える
np.random.seed(42)  # 再現性のため
noise = np.random.normal(0, 5, size=extended_population.shape)  # 平均0、標準偏差5のノイズ
noisy_population = extended_population + noise

# 前半部分と後半部分を分ける
half_index = len(extended_years) // 2
train_years = extended_years[:half_index]  # 前半の年度データ
train_population = noisy_population[:half_index]  # ノイズを加えた前半の人口データ

# 差分近似で微分値を計算する関数（新しい近似式）
def calculate_derivatives_new(years, population):
    dN_dt = []
    for i in range(1, len(population) - 1):
        # 新しい近似式を適用
        derivative = 0.5 * (
            (population[i + 1] - population[i]) / (years[i + 1] - years[i]) + 
            (population[i] - population[i - 1]) / (years[i] - years[i - 1])
        )
        dN_dt.append(derivative)
    return np.array(dN_dt)

# 新しい微分値と中間のNデータ
dN_dt_train_new = calculate_derivatives_new(train_years, train_population)
N_mid_train_new = train_population[1:-1]  # 中間の人口値

# 損失関数の定義（新しい近似式を使用）
def loss_function_new(params):
    r, C = params
    model_derivative = r * N_mid_train_new * (1 - N_mid_train_new / C)
    return np.sum((dN_dt_train_new - model_derivative) ** 2)

# 初期値と境界条件
initial_guess = [0.02, 200]  # rとCの初期値
bounds = [(0, None), (0, None)]  # rとCは正の値のみ

# 最適化
result_new = minimize(loss_function_new, initial_guess, bounds=bounds)

# 最適化結果
r_opt_new, C_opt_new = result_new.x
min_loss_new = result_new.fun  # 最小化された損失関数の値
print(f"Estimated r (new): {r_opt_new}")
print(f"Estimated C (new): {C_opt_new}")
print(f"Minimized loss (new): {min_loss_new}")

# 推定されたモデルを差分方程式でシミュレーション
def simulate_population_difference(t, N0, r, C):
    N_sim = [N0]
    dt = t[1] - t[0]
    for i in range(len(t) - 1):
        N_current = N_sim[-1]
        dN = r * N_current * (1 - N_current / C) * dt
        N_sim.append(N_current + dN)
    return np.array(N_sim)

# 前半のデータを使用してシミュレーション
predicted_population = simulate_population_difference(extended_years, noisy_population[0], r_opt_new, C_opt_new)

# グラフの表示
plt.figure(figsize=(10, 6))

# ノイズを加えたデータ（前半と後半に分けて表示）
plt.plot(years, population, 'o', label="Original Data (Observed)", color='blue')
plt.plot(extended_years, noisy_population, '--', label="Spline Interpolation + Noise (Extended)", color='orange')

# シミュレーション結果を小さな丸で表示
plt.scatter(extended_years, predicted_population, color='green', label="Simulation (Predicted)", s=10)

# 境目を示す縦線（前半と後半の境界）
plt.axvline(x=extended_years[half_index], color='red', linestyle='--', label="Boundary between Train and Test Data")

# ラベルとタイトルの設定
plt.xlabel("Year")
plt.ylabel("Population (×10^6)")
plt.title("Population Data: Observed, Extended with Noise, and Simulation Comparison")
plt.legend()
plt.grid()

# 見やすいようにマーカーを追加（元データとシミュレーションデータを強調）
plt.scatter(years, population, color='blue', zorder=5)
plt.scatter(extended_years, predicted_population, color='green', zorder=5, s=10)

plt.show()

グラフはこんな形でホームランしてしまった。どうしてだろうか。

耐ノイズへの変形

# 差分近似で微分値を計算する関数（新しい近似式）
def calculate_derivatives_new(years, population):
    dN_dt = []
    for i in range(1, len(population) - 1):
        # 新しい近似式を適用
        derivative = 0.5 * (
            (population[i + 1] - population[i]) / (years[i + 1] - years[i]) + 
            (population[i] - population[i - 1]) / (years[i] - years[i - 1])
        )
        dN_dt.append(derivative)
    return np.array(dN_dt)

# 新しい微分値と中間のNデータ
dN_dt_train_new = calculate_derivatives_new(train_years, train_population)
N_mid_train_new = train_population[1:-1]  # 中間の人口値

# 損失関数の定義（新しい近似式を使用）
def loss_function_new(params):
    r, C = params
    model_derivative = r * N_mid_train_new * (1 - N_mid_train_new / C)
    return np.sum((dN_dt_train_new - model_derivative) ** 2)

# 初期値と境界条件
initial_guess = [0.02, 200]  # rとCの初期値
bounds = [(0, None), (0, None)]  # rとCは正の値のみ

# 最適化
result_new = minimize(loss_function_new, initial_guess, bounds=bounds)

コードをよく見ると気になることがあった。minimizeに渡している関数はloss_function_newという関数である。この関数ではdN_dt_train_newという平均微分の値を格納しており、微分の値が厳密に一致するように最適化している。ここに問題がある。
　今、一番最適化して予測したいのは人口であり、微分の値ではない。プロンプトの出力結果を見ると、本当はちゃんと書いてあった。ChatGPT先生は損失関数を微分の二乗誤差として定義していた。

なので、プロンプトを流し込む

損失関数の定義を変更してください。
今は、損失関数はdN/dtに対する誤差として定義していますが、Nに対する誤差として定義してください。

耐ノイズへの変形のシミュレーション

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.interpolate import CubicSpline

# 年度と人口データ（元データ）
years = np.array([1820, 1830, 1840, 1850, 1860, 1870, 1880, 1890, 1900, 1910, 1920, 1930])
population = np.array([9.6, 12.9, 17.1, 23.2, 31.4, 38.6, 50.2, 62.9, 76.0, 92.0, 106.5, 123.2])

# 3次スプライン補完を行う
spline = CubicSpline(years, population)

# 新しい1年刻みの年度データを作成
extended_years = np.arange(years[0], years[-1] + 1, 1)

# 1年刻みで補完した人口データを計算
extended_population = spline(extended_years)

# N(0,5)の白色ノイズを加える
np.random.seed(42)  # 再現性のため
noise = np.random.normal(0, 5, size=extended_population.shape)  # 平均0、標準偏差5のノイズ
noisy_population = extended_population + noise

# 前半部分と後半部分を分ける
half_index = len(extended_years) // 2
train_years = extended_years[:half_index]  # 前半の年度データ
train_population = noisy_population[:half_index]  # ノイズを加えた前半の人口データ

# シミュレーション関数（人口データから予測人口を計算）
def simulate_population_difference(t, N0, r, C):
    N_sim = [N0]
    dt = t[1] - t[0]
    for i in range(len(t) - 1):
        N_current = N_sim[-1]
        dN = r * N_current * (1 - N_current / C) * dt
        N_sim.append(N_current + dN)
    return np.array(N_sim)

# 損失関数の定義（Nに対する誤差として定義）
def loss_function_new(params):
    r, C = params
    # シミュレーションを行い、前半部分の人口を計算
    predicted_population = simulate_population_difference(extended_years, noisy_population[0], r, C)
    
    # 前半部分の予測人口と実際の人口との誤差（二乗和）を計算
    error = np.sum((predicted_population[:half_index] - noisy_population[:half_index]) ** 2)
    return error

# 最適化中にパラメーターと誤差を表示するためのコールバック関数
param_history = []
loss_history = []

def track_progress(params):
    r, C = params
    loss = loss_function_new(params)
    param_history.append((r, C))
    loss_history.append(loss)
    print(f"r: {r:.5f}, C: {C:.5f}, Loss: {loss:.5f}")

# 初期値と境界条件
initial_guess = [0.02, 200]  # rとCの初期値
bounds = [(0, None), (0, None)]  # rとCは正の値のみ

# 最適化
result_new = minimize(loss_function_new, initial_guess, bounds=bounds, callback=track_progress)

# 最適化結果
r_opt_new, C_opt_new = result_new.x
min_loss_new = result_new.fun  # 最小化された損失関数の値
print(f"Estimated r (new): {r_opt_new}")
print(f"Estimated C (new): {C_opt_new}")
print(f"Minimized loss (new): {min_loss_new}")

# シミュレーション結果を再計算
predicted_population = simulate_population_difference(extended_years, noisy_population[0], r_opt_new, C_opt_new)

# 最適化過程のパラメータと損失の変化をプロット
plt.figure(figsize=(10, 6))

# パラメータの履歴を表示
param_history = np.array(param_history)
plt.subplot(2, 1, 1)
plt.plot(param_history[:, 0], label="r (growth rate)")
plt.plot(param_history[:, 1], label="C (carrying capacity)")
plt.title("Optimization Progress: Parameter History")
plt.xlabel("Iteration")
plt.ylabel("Parameter Value")
plt.legend()

# 損失の履歴を表示
plt.subplot(2, 1, 2)
plt.plot(loss_history, label="Loss (Objective Function)")
plt.title("Optimization Progress: Loss History")
plt.xlabel("Iteration")
plt.ylabel("Loss Value")
plt.legend()

plt.tight_layout()
plt.show()

# グラフの表示
plt.figure(figsize=(10, 6))

# ノイズを加えたデータ（前半と後半に分けて表示）
plt.plot(years, population, 'o', label="Original Data (Observed)", color='blue')
plt.plot(extended_years, noisy_population, '--', label="Spline Interpolation + Noise (Extended)", color='orange')

# シミュレーション結果を小さな丸で表示
plt.scatter(extended_years, predicted_population, color='green', label="Simulation (Predicted)", s=10)

# 境目を示す縦線（前半と後半の境界）
plt.axvline(x=extended_years[half_index], color='red', linestyle='--', label="Boundary between Train and Test Data")

# ラベルとタイトルの設定
plt.xlabel("Year")
plt.ylabel("Population (×10^6)")
plt.title("Population Data: Observed, Extended with Noise, and Simulation Comparison")
plt.legend()
plt.grid()

# 見やすいようにマーカーを追加（元データとシミュレーションデータを強調）
plt.scatter(years, population, color='blue', zorder=5)
plt.scatter(extended_years, predicted_population, color='green', zorder=5, s=10)

plt.show()

r: 0.02400, C: 200.00000, Loss: 1463.90934
r: 0.02534, C: 200.00000, Loss: 1427.89472
r: 0.02505, C: 200.00000, Loss: 1425.14808
r: 0.02506, C: 200.00000, Loss: 1425.14211
r: 0.02506, C: 200.00000, Loss: 1425.14211
Estimated r (new): 0.02506146162474434
Estimated C (new): 200.0000000790879
Minimized loss (new): 1425.1421080398673

ぼちぼちうまいことフィッティングできているようだ。

差分方程式によるシミュレーション

先ほど耐ノイズへの変形を行ったが、そのコードを見て思ったことがある。

# シミュレーション関数（人口データから予測人口を計算）
def simulate_population_difference(t, N0, r, C):
    N_sim = [N0]
    dt = t[1] - t[0]
    for i in range(len(t) - 1):
        N_current = N_sim[-1]
        dN = r * N_current * (1 - N_current / C) * dt
        N_sim.append(N_current + dN)
    return np.array(N_sim)

# 損失関数の定義（Nに対する誤差として定義）
def loss_function_new(params):
    r, C = params
    # シミュレーションを行い、前半部分の人口を計算
    predicted_population = simulate_population_difference(extended_years, noisy_population[0], r, C)
    
    # 前半部分の予測人口と実際の人口との誤差（二乗和）を計算
    error = np.sum((predicted_population[:half_index] - noisy_population[:half_index]) ** 2)
    return error

微分の導出がなくなっている。確かにデータとして不要になってしまった。そこで、はたと思いついたのが差分方程式である。

このように書くことで二項間漸化式に帰着できる。これを用いて、シミュレーションしてみよう。この場合、r,C,N[0],N[1]を最適化パラメータとした。

差分方程式によるシミュレーションのコード

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.interpolate import CubicSpline

# 年度と人口データ（元データ）
years = np.array([1820, 1830, 1840, 1850, 1860, 1870, 1880, 1890, 1900, 1910, 1920, 1930])
population = np.array([9.6, 12.9, 17.1, 23.2, 31.4, 38.6, 50.2, 62.9, 76.0, 92.0, 106.5, 123.2])

# 3次スプライン補完を行う
spline = CubicSpline(years, population)

# 新しい1年刻みの年度データを作成
extended_years = np.arange(years[0], years[-1] + 1, 1)

# 1年刻みで補完した人口データを計算
extended_population = spline(extended_years)

# N(0,5)の白色ノイズを加える
np.random.seed(42)  # 再現性のため
noise = np.random.normal(0, 5, size=extended_population.shape)  # 平均0、標準偏差5のノイズ
noisy_population = extended_population + noise

# 前半部分と後半部分を分ける
half_index = len(extended_years) // 2
train_years = extended_years[:half_index]  # 前半の年度データ
train_population = noisy_population[:half_index]  # ノイズを加えた前半の人口データ

# 差分方程式のモデル
def diff_eq_model(params, t):
    r, C, N0, N1 = params  # r, C, N0, N1をパラメータとして受け取る
    dt = np.diff(t)  # 時間の間隔
    N = np.zeros(len(t))
    N[0] = N0
    N[1] = N1  # N1を設定
    for i in range(1, len(t) - 1):
        N[i+1] = N[i-1] + dt[i-1] * r * N[i] * (1 - N[i] / C)
    return N

# 目的関数（最小二乗誤差）
def objective(params, t, population_data):
    r, C, N0, N1 = params
    if N0 < 0 or N1 < 0:  # N0とN1は0以上でなければならないので、制約を追加
        return np.inf
    predicted_population = diff_eq_model(params, t)
    error = np.sum((predicted_population[1:] - population_data[1:])**2)
    return error

# 最適化中にパラメーターと誤差を表示するためのコールバック関数
param_history = []
loss_history = []

def track_progress(params):
    r, C, N0, N1 = params
    loss = objective(params, train_years, train_population)
    param_history.append((r, C, N0, N1))
    loss_history.append(loss)
    print(f"r: {r:.5f}, C: {C:.5f}, N0: {N0:.5f}, N1: {N1:.5f}, Loss: {loss:.5f}")

# 最適化の実行
initial_guess = [0.1, 50, train_population[0], train_population[1]]  # r=0.1, C=50, 初期人口N0, N1
result = minimize(objective, initial_guess, args=(train_years, train_population), bounds=((0, None), (0, None), (0, None), (0, None)), callback=track_progress)

# 最適化されたr, C, N0, N1を取得
r_opt, C_opt, N0_opt, N1_opt = result.x
print(f"最適化されたr: {r_opt}")
print(f"最適化されたC: {C_opt}")
print(f"最適化されたN0: {N0_opt}")
print(f"最適化されたN1: {N1_opt}")

# 最適化されたパラメータで予測を実行
predicted_population = diff_eq_model([r_opt, C_opt, N0_opt, N1_opt], extended_years)

# 最適化過程のパラメータと損失の変化をプロット
param_history = np.array(param_history)
plt.figure(figsize=(10, 6))

# パラメータの履歴を表示
plt.subplot(2, 1, 1)
plt.plot(param_history[:, 0], label="r (growth rate)")
plt.plot(param_history[:, 1], label="C (carrying capacity)")
plt.plot(param_history[:, 2], label="N0 (initial population)")
plt.plot(param_history[:, 3], label="N1 (1st population)")
plt.title("Optimization Progress: Parameter History")
plt.xlabel("Iteration")
plt.ylabel("Parameter Value")
plt.legend()

# 損失の履歴を表示
plt.subplot(2, 1, 2)
plt.plot(loss_history, label="Loss (Objective Function)")
plt.title("Optimization Progress: Loss History")
plt.xlabel("Iteration")
plt.ylabel("Loss Value")
plt.legend()

plt.tight_layout()
plt.show()

# グラフの表示
plt.figure(figsize=(10, 6))

# ノイズを加えたデータ（前半と後半に分けて表示）
plt.plot(years, population, 'o', label="Original Data (Observed)", color='blue')
plt.plot(extended_years, noisy_population, '--', label="Spline Interpolation + Noise (Extended)", color='orange')

# シミュレーション結果を小さな丸で表示
plt.scatter(extended_years, predicted_population, color='green', label="Simulation (Predicted)", s=10)

# 境目を示す縦線（前半と後半の境界）
plt.axvline(x=extended_years[half_index], color='red', linestyle='--', label="Boundary between Train and Test Data")

# ラベルとタイトルの設定
plt.xlabel("Year")
plt.ylabel("Population (×10^6)")
plt.title("Population Data: Observed, Extended with Noise, and Simulation Comparison")
plt.legend()
plt.grid()

# 見やすいようにマーカーを追加（元データとシミュレーションデータを強調）
plt.scatter(years, population, color='blue', zorder=5)
plt.scatter(extended_years, predicted_population, color='green', zorder=5, s=10)

# グラフを表示
plt.show()

(略)
r: 0.05783, C: 8280.50355, N0: 9.27658, N1: 9.27114, Loss: 1103.87479
r: 0.05783, C: 8280.50355, N0: 9.27662, N1: 9.27121, Loss: 1103.87450
r: 0.05783, C: 8280.50355, N0: 9.27650, N1: 9.27100, Loss: 1103.87443
r: 0.05784, C: 8280.50355, N0: 9.27551, N1: 9.26927, Loss: 1103.87402
r: 0.05786, C: 8280.50355, N0: 9.27388, N1: 9.26643, Loss: 1103.87353
r: 0.05786, C: 8280.50355, N0: 9.27240, N1: 9.26385, Loss: 1103.87321
r: 0.05787, C: 8280.50355, N0: 9.27204, N1: 9.26322, Loss: 1103.87314
r: 0.05786, C: 8280.50355, N0: 9.27212, N1: 9.26335, Loss: 1103.87313
r: 0.05786, C: 8280.50355, N0: 9.27215, N1: 9.26341, Loss: 1103.87313
最適化されたr: 0.05786402532275687
最適化されたC: 8280.503551306912
最適化されたN0: 9.272148323918822
最適化されたN1: 9.263410983850378

学習データにはきれいにフィットしたが、長期予測がオーバーシュートしている傾向にある。

感想

「もうちょいきれいにフィッティングできねぇかなぁ。」

ノイズの無いデータに対してでも割と長期予測は外れてるので、もうちょい頑張れないかなぁ。感はありました。一方で

「実験設定失敗したなぁ。」

という感触もあります。なんとなく本に書いてあるデータを参考にすれば、ばっちりフィッティングできるだろうというヤマカンで始めましたが、よくよく考えてみるとそうではなかった。もともとの人口増加のデータがヴェアフルストの人口モデルの微分方程式にフィッティングできるとは限らなかった。そう考えると、まず初めにヴェアフルストの人口モデルでデータを生成したのち、そのデータに対してフィッティングして、精度をどうのこうの議論するのが筋だったな。というのは思うことでした。

「シミュレーションは昔勉強したなぁ」

大学時代数値計算の授業があり、その時に、オイラー法なり、ルンゲクッタ法なりを実装した記憶があり、どこで使うねん。と思いました。今使ってる。そう。この問題って、実は微分方程式のシミュレーションにもかかわる話なので、ここも結構気にしないといけない部分で、陰解法・陽解法いろいろあったな・・・みたいなことを思い出しました。

「ChatGPTすごく便利」

ChatGPTかなり書けるな。というイメージはありました。だから、ほとんどChatGPTと対話するだけで検証をしてみよう。という試みでした。今回、エディターは1つも使ってなくて、ChatGPTとPython実行用の環境としてGoogleColabを使っただけです。これだけでここまで検証できた。ChatGPTの便利な使い方として、「グラフを書く」というのが個人的に好きな使い方です。毎回毎回可視化のためにmatplotlibのドキュメントを見に行くのが死ぬほどつらかったので、

最適化した結果の誤差を評価したいです。
元のデータと、推定された方程式のグラフを表示して欲しいです。

ぐらいのノリでグラフを描画するコードを出してくれるのは、超楽です。一方で、問題だなぁと思うのが、議論が脇道にそれた後にコーディングに戻ってきて欲しいときとか、一定時間たった後にChatGPTに問い合わせたりとか、そういう場合に、元のコードや議論を忘れていたりすることがあって、結構厄介でした。

「もうちょい楽にならんかな」

という気持ちも少しあります。これは、ヴェアフルストの人口モデルに関して、えっさほいさとやっていましたが、任意の微分方程式に対して、こういう試行ができるフレームワークあると便利なんだけど、むずいだろうなぁと思いました。

「最適化できてなくない？」

みたいな感触もあります。もうちょいパラメータフィッティングがんばれるんじゃない？感はあります。関数の形などは利用せずにフィッティングしているので、その意味ではまだ頑張れるんじゃない？感はあります。あとAdamとかああいうモダンなやつってscipyで使えないのかなぁ。とかぼんやりと思いました。