既知モデル＋未知項（Grey-box）で「モデル破綻の理由」を掘る
研究でよくある状況：
- 既知の理論モデル（微分方程式・反応式・経験式）がある
- でも、データに当てると途中から外れる
- 「モデルがダメ」なのか「測定がダメ」なのか「前提が破れている」のか分からない
このとき、いきなりブラックボックス（巨大NN）に逃げずに、
既知モデル + 未知項（Grey-box / Hybrid / Model-discrepancy）
として「どこが足りないのか」を“形”として掘ると、次の一手（仮説・追加実験・モデル改良）が打ちやすくなります。
この記事では、初心者向けに ダミーデータで
- 既知モデル（指数成長）を当てる → 破綻する
- 残差（未知項）を学習して 破綻の理由を可視化
- Grey-box でモデルを修正 → フィットが改善
までを Google Colabで動くコードで再現します。
TL;DR
- Grey-box = 「信じる部分（既知モデル）」は残して「足りない部分（未知項）」だけ学習
- “モデル破綻”は、未知項の 時間依存 / 状態依存 / 閾値依存 の形として現れやすい
- まずは未知項を 単純な形（多項式・周期・区分）で当てて、原因仮説を作るのが強い
（必要なら後でNN等に拡張）
1. 用語整理
1.1 Grey-box（グレーボックス）とは？
「白（ホワイト）＝理論モデル100%」と「黒（ブラック）＝データ駆動100%」の中間です。
- 既知モデル：理論・物理・化学・生理・工学の理解を反映
- 未知項：既知モデルが説明しきれない“ずれ”を吸収（ただし、吸収した形を解析して理由を掘る）
典型形はこうです：
- 既知：
du/dt = f_known(u, t; θ) - Grey-box：
du/dt = f_known(u, t; θ) + g_unknown(u, t; φ)
g_unknown を「残差（model discrepancy）」として学習し、その形を見ることで、
「何が足りないのか」を推理します。
1.2 「モデル破綻の理由」を掘るとは？
観測 u(t) が既知モデルから外れるとき、原因候補はだいたい次のどれかです：
-
未モデリングの入力（外力・温度・周期・操作量）：
g(t)が周期的・ステップ的になる -
非線形（飽和・閾値・相転移・阻害）：
g(u)がu^2やS字・閾値依存になる - 観測モデルの破綻（センサー飽和、遅れ、オフセット）：残差が観測値と相関したり、特定レンジで急変する
- 条件の混在（バッチ差・個体差）：群で残差が分かれる
今回は「指数成長モデルが破綻する理由」として、飽和（資源制限）を例にします。
2. 今回のトイ問題：指数成長モデルが途中から外れる理由は「飽和」という例
よくある勘違いの例：
-
信じたいモデル（既知）：指数成長
du/dt = r u -
実際（真の生成）：ロジスティック成長（飽和）
du/dt = r u (1 - u/K) = r u - (r/K) u^2
指数モデルでフィットすると、初期は合うのに、途中から
- 予測が 大きくなりすぎる
- 残差が 系統的に出る
という典型パターンになります。
この「系統的残差」の形を見て、
未知項が -const * u^2 っぽい → 「飽和が抜けてる」
→ じゃあロジスティックに改良しよう
という仮説を作るのが本記事のゴールです。
3. Google Colabで実行可能なPythonコード
3.1 実行の流れ
- Colabを開く（CPUでも可）
- 以下のセルを上から順に実行
-
fig/に画像が保存される
セル1：準備＆ダミーデータ作成（真の生成はロジスティック）
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
from sklearn.linear_model import LinearRegression, Ridge
os.makedirs("fig", exist_ok=True)
np.random.seed(0)
# --- True dynamics (toy) ---
r_true = 1.2
K_true = 1.0
u0_true = 0.05
t_end = 10.0
n = 201
t = np.linspace(0, t_end, n)
dt = t[1] - t[0]
def logistic_solution(t, r, K, u0):
# analytic solution of logistic equation
return K / (1.0 + (K/u0 - 1.0) * np.exp(-r*t))
u_true = logistic_solution(t, r_true, K_true, u0_true)
# noisy observations
sigma = 0.03
u_obs = u_true + np.random.randn(n) * sigma
# clip for safety (log / smoothing)
u_obs_clip = np.clip(u_obs, 1e-4, None)
plt.figure(figsize=(7.2,4.2))
plt.plot(t, u_true, label="ground truth (logistic)")
plt.scatter(t, u_obs, s=12, alpha=0.7, label="noisy obs")
plt.xlabel("t")
plt.ylabel("u(t)")
plt.title("Toy data: true = logistic (saturation), but we will fit exponential first")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig1_toy_data.png", dpi=160)
plt.show()
print("Saved: fig/fig1_toy_data.png")
セル2：既知モデル（指数成長）をフィットして「破綻」を確認
指数成長 u(t)=u0*exp(r t) なので、初期の範囲では
log u ≈ log u0 + r t
として線形回帰できます。
# Fit exponential model on early time window (where u is small)
t_fit_max = 2.0
mask_early = (t <= t_fit_max)
X = t[mask_early].reshape(-1, 1)
y = np.log(u_obs_clip[mask_early])
lin = LinearRegression().fit(X, y)
r_exp_hat = float(lin.coef_[0])
log_u0_hat = float(lin.intercept_)
u0_hat = float(np.exp(log_u0_hat))
print(f"Estimated exponential params (early fit): r_hat={r_exp_hat:.3f}, u0_hat={u0_hat:.3f}")
u_exp_pred = u0_hat * np.exp(r_exp_hat * t)
plt.figure(figsize=(7.2,4.2))
plt.plot(t, u_true, label="ground truth (logistic)")
plt.scatter(t, u_obs, s=12, alpha=0.7, label="noisy obs")
plt.plot(t, u_exp_pred, linewidth=2.5, label="known model: exponential fit (early only)")
plt.xlabel("t")
plt.ylabel("u(t)")
plt.title("Known model breaks: exponential overshoots after saturation")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig2_exponential_breakdown.png", dpi=160)
plt.show()
print("Saved: fig/fig2_exponential_breakdown.png")
読み方（重要）
- 初期（t<=2）は合う
- でも後半で指数モデルが 上に突き抜ける
→ 「指数成長だけでは説明できない」サイン
セル3：残差（未知項）を作る：まずは du/dt を推定して “微分の残差” を見る
Grey-boxの基本は
du/dt ≈ f_known(u,t) + g_unknown(u,t)
なので、データから du/dt を推定し、既知項を引いて残差を作ります。
注意：微分はノイズに弱いので、ここでは Savitzky-Golay で平滑化＋微分をします。
実データでは、窓幅や手法（スプライン/GP/カルマン）を工夫してください。
# choose robust window length (odd)
win = 31
if win >= n:
win = n-1 if (n-1) % 2 == 1 else n-2
if win < 7:
win = 7
poly = 3
u_smooth = savgol_filter(u_obs_clip, window_length=win, polyorder=poly)
du_dt = savgol_filter(u_obs_clip, window_length=win, polyorder=poly, deriv=1, delta=dt)
# known derivative under exponential model: f_known = r_hat * u
f_known = r_exp_hat * u_smooth
# discrepancy (unknown term proxy): g_hat_raw(t) = du/dt - f_known
g_raw = du_dt - f_known
plt.figure(figsize=(7.2,4.2))
plt.plot(t, g_raw, label="raw discrepancy g(t) = du/dt - r_hat*u")
plt.axhline(0, linestyle="--")
plt.xlabel("t")
plt.ylabel("discrepancy (approx)")
plt.title("Discrepancy over time (hint: systematic, not random)")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig3_discrepancy_over_time.png", dpi=160)
plt.show()
print("Saved: fig/fig3_discrepancy_over_time.png")
セル4：残差の“形”を見る：g は u に依存してない？（飽和なら u^2 が出る）
モデル破綻の理由を掘るコツは、
- 残差を 時間 t に対して見る（外力・周期・操作）
- 残差を 状態 u に対して見る（非線形・飽和・閾値）
- それでもダメなら 外部特徴量 x を入れる（温度・濃度・装置ログ等）
今回は飽和なので g(u) を見ます。
plt.figure(figsize=(6.4,4.8))
plt.scatter(u_smooth, g_raw, s=12, alpha=0.7)
plt.axhline(0, linestyle="--")
plt.xlabel("u (smoothed)")
plt.ylabel("g_raw = du/dt - r_hat*u")
plt.title("Discrepancy vs state u (shape reveals missing physics)")
plt.tight_layout()
plt.savefig("fig/fig4_discrepancy_vs_u.png", dpi=160)
plt.show()
print("Saved: fig/fig4_discrepancy_vs_u.png")
期待される読み（今回のトイでは）
- u が大きくなるほど g が 負に寄る
- しかも曲線っぽい → u^2 っぽい非線形が怪しい
セル5：未知項を “まずは単純に” 当てる（多項式：a0 + a1 u + a2 u^2）
ここが重要ポイント：
いきなりNNにせず、まずは 解釈しやすい形で当てます。
g(u) ≈ a0 + a1 u + a2 u^2
そして
-
du/dt = r_hat u + g(u)
= (r_hat + a1) u + a2 u^2 + a0
となるので、
r_grey = r_hat + a1-
a2 < 0なら飽和っぽい - ロジスティックに近いなら
a2 ≈ -(r/K)→K ≈ -r_grey/a2
が推測できます（※トイだからできる“推理”）。
# Fit g(u) with a simple polynomial (Ridge for stability)
U = u_smooth.reshape(-1, 1)
X_poly = np.hstack([np.ones_like(U), U, U**2]) # [1, u, u^2]
ridge = Ridge(alpha=1e-3, fit_intercept=False).fit(X_poly, g_raw)
a0, a1, a2 = ridge.coef_
r_grey = r_exp_hat + a1
print("Estimated discrepancy model: g(u) = a0 + a1*u + a2*u^2")
print(f"a0={a0:.4f}, a1={a1:.4f}, a2={a2:.4f}")
print(f"-> implied r_grey = r_hat + a1 = {r_grey:.4f}")
K_hat = None
if a2 < 0:
K_hat = - r_grey / a2
print(f"-> implied K_hat ≈ -r_grey/a2 = {K_hat:.3f} (toy true K={K_true})")
else:
print("a2 is not negative -> saturation is not supported by this simple fit (check data / window / model).")
# visualize fitted g(u)
u_line = np.linspace(u_smooth.min(), u_smooth.max(), 200)
g_fit = a0 + a1*u_line + a2*(u_line**2)
plt.figure(figsize=(6.4,4.8))
plt.scatter(u_smooth, g_raw, s=10, alpha=0.5, label="raw discrepancy")
plt.plot(u_line, g_fit, linewidth=3, label="fitted g(u) (quadratic)")
plt.axhline(0, linestyle="--")
plt.xlabel("u")
plt.ylabel("g(u)")
plt.title("Learned unknown term g(u): interpretable shape")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig5_g_of_u.png", dpi=160)
plt.show()
print("Saved: fig/fig5_g_of_u.png")
ここまでで「モデル破綻の理由」にかなり踏み込めます：
-
a2が負で、g(u)が-u^2的
→ 高濃度（大きいu）ほど成長が抑制される（飽和・阻害）
→ 指数成長の前提（無限の資源）が破れている
セル6：Grey-boxモデルで再シミュレーションして「破綻が直るか」確認
Grey-box ODE を
du/dt = r_hat*u + (a0 + a1*u + a2*u^2)
として数値積分（RK4）します。
def rk4(f, y0, t):
y = np.zeros_like(t, dtype=float)
y[0] = float(y0)
for i in range(len(t)-1):
dt = t[i+1] - t[i]
ti = t[i]
yi = y[i]
k1 = f(yi, ti)
k2 = f(yi + 0.5*dt*k1, ti + 0.5*dt)
k3 = f(yi + 0.5*dt*k2, ti + 0.5*dt)
k4 = f(yi + dt*k3, ti + dt)
y[i+1] = yi + dt*(k1 + 2*k2 + 2*k3 + k4)/6.0
return y
def f_grey(u, t_):
return r_exp_hat*u + (a0 + a1*u + a2*u*u)
u0_sim = float(u_obs_clip[0]) # use first obs as initial value
u_grey_pred = rk4(f_grey, u0_sim, t)
plt.figure(figsize=(7.2,4.2))
plt.plot(t, u_true, label="ground truth (logistic)")
plt.scatter(t, u_obs, s=12, alpha=0.7, label="noisy obs")
plt.plot(t, u_exp_pred, linewidth=2.2, label="known exponential (breaks)")
plt.plot(t, u_grey_pred, linewidth=3.0, label="grey-box: exponential + learned g(u)")
plt.xlabel("t")
plt.ylabel("u(t)")
plt.title("Grey-box repairs breakdown and hints the missing term")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig6_greybox_repair.png", dpi=160)
plt.show()
print("Saved: fig/fig6_greybox_repair.png")
4. 結果の解釈：この例で「モデル破綻の理由」は何だったのか？
今回のトイでは、指数モデルが破綻した理由はシンプルで：
- 本当は
du/dt = r u - (r/K) u^2なのに -
du/dt = r uしか入れていなかった
つまり、“飽和（資源制限）”という非線形項が抜けていたのが原因です。
Grey-boxで g(u) を学習すると：
-
g(u)が-u^2的に下がる - 「大きい u ほど成長が抑制される」という構造が可視化される
→ そこから
- ロジスティック
- 阻害項（Hill型）
- 枯渇モデル
など「次の候補モデル」を提案できます。
5. 実データに適用するときの“型”
5.1 手順（おすすめ）
- 既知モデルを当てる（当たる範囲・外れる範囲を確認）
- データから
du/dt（や差分）を推定し、残差g = du/dt - f_knownを作る - 残差をまず 可視化する
-
g(t)：周期・ステップ・ドリフトがあるか -
g(u)：非線形・閾値があるか -
g(x)：温度・操作・装置ログと相関があるか
-
- いきなりNNにせず、まずは 単純な形で当てる（解釈のため）
- 多項式、周期（sin/cos）、区分線形、スプライン、L1で疎に（SINDy風）
- 仮説が立ったら
- 追加実験（どこを測るべきか）
- モデル改良（項を足す、パラメータ推定し直す）
を回す
5.2 よくある落とし穴
- 微分がノイズで壊れる：平滑化・スプライン・GP・カルマンを検討
- 未知項がノイズを吸いすぎる：正則化・単純モデル優先・検証用データを分ける
- “当たった形”の過解釈：外部変数や実験条件と突き合わせる（因果に注意）
6. もう一歩：この考え方の発展形（興味があれば）
- Known + NN（Neural ODE / UDE）：未知項をNNにして、微分方程式を微分可能に学習
- SINDy（スパース回帰）：未知項を「少数の項」で説明する（解釈強い）
- PINN系：観測が少ない／PDE系で、物理制約を強くかけながら未知項も推定
- 実験計画（次に測る点）：未知項が大きい領域・感度が高い領域を狙って追加測定
ただし、まずは本記事のように
“残差の形を見て、単純に当てて、原因仮説を作る”
が汎用的初期検討としておすすめです。
まとめ
- Grey-box（既知モデル＋未知項）は、モデルを捨てずに 破綻理由を掘るための武器
- 残差
gを t/u/外部変数に対して可視化すると、欠落要因の仮説が立つ - 今回は
g(u) ~ -u^2→ 飽和が抜けていたと見抜ける例を示した - 研究では「モデル改良」と「追加測定」の意思決定に直結する