はじめに
「指数関数的減衰」という言葉は、物理、化学、生物、医学、工学、データ解析など、かなり広い分野で出てきます。
たとえば、次のような現象です。
- 温度差が時間とともに小さくなる
- 蛍光強度が時間とともに弱くなる
- 薬物濃度が体内で減っていく
- コンデンサの電圧が放電で下がる
- ある応答が刺激後にゆっくり元へ戻る
これらはすべて、状況によっては「指数関数的減衰」で近似できます。
この記事では、指数関数的減衰について、用語の説明から始めて、
- 時定数とは何か
- 半減期とは何か
- 時定数と半減期はどう関係するのか
- ノイズ入りのダミーデータから指数関数的減衰をフィッティングするにはどうするのか
を、Pythonで図を作りながら確認します。
実データではなく、まずは仕組みを理解するためのダミーデータを使います。
この記事のゴール
この記事のゴールは、次の式を「ただの式」ではなく、「図とコードで意味が分かる式」として理解することです。
$$
y(t) = C + A e^{-t/\tau}
$$
ここで、
| 記号 | 意味 |
|---|---|
| $t$ | 時間 |
| $y(t)$ | 時刻 $t$ における測定値 |
| $A$ | 初期振幅。ベースラインからどれくらい上にあるか |
| $C$ | ベースライン。十分時間が経った後に近づく値 |
| $\tau$ | 時定数。減衰の速さを決める値 |
| $e$ | ネイピア数。約 $2.718$ |
ポイントは、$y(t)$ そのものがゼロに向かうとは限らないことです。
多くの実験データでは、信号はゼロではなく、あるベースライン $C$ に近づいていきます。そのため、
$$
y(t) - C = A e^{-t/\tau}
$$
と考えると分かりやすくなります。
指数関数的減衰とは何か
指数関数的減衰は、ざっくり言うと、
残っている量に比例して減っていく現象
です。
「毎秒1だけ減る」のではありません。
そうではなく、
残っている量のうち、一定の割合が減っていく
というイメージです。
たとえば、最初に100あったものが、一定時間ごとに半分になるなら、
100 → 50 → 25 → 12.5 → 6.25 → ...
のように減ります。
減る量そのものは小さくなっていきますが、「割合」としては同じように減っていきます。
これが指数関数的減衰の直感です。
時定数とは何か
指数関数的減衰でよく出てくるのが、時定数 $\tau$ です。
式は次の通りです。
$$
y(t) - C = A e^{-t/\tau}
$$
ここで $t=\tau$ を代入すると、
$$
y(\tau) - C = A e^{-1}
$$
となります。
$e^{-1}$ は約 $0.368$ なので、
$$
y(\tau)-C \approx 0.368 A
$$
です。
つまり、時定数 $\tau$ は、
ベースラインからの差が、最初の約36.8%まで減る時間
です。
ここで注意したいのは、時定数は「半分になる時間」ではない、ということです。
半減期とは何か
半減期 $t_{1/2}$ は、その名の通り、
ベースラインからの差が半分になる時間
です。
つまり、
$$
y(t_{1/2}) - C = \frac{A}{2}
$$
となる時間です。
指数関数的減衰の式に代入すると、
$$
\frac{A}{2} = A e^{-t_{1/2}/\tau}
$$
両辺を $A$ で割ると、
$$
\frac{1}{2} = e^{-t_{1/2}/\tau}
$$
両辺の自然対数を取ると、
$$
\log \frac{1}{2} = -\frac{t_{1/2}}{\tau}
$$
したがって、
$$
t_{1/2} = \tau \log 2
$$
となります。
つまり、半減期と時定数の関係は、
$$
t_{1/2} \approx 0.693 \tau
$$
です。
半減期は時定数より少し短い時間になります。
Google Colabで実行するコード
以下のコードをGoogle Colabで実行してください。
このコードでは、
- 時定数と半減期の関係を図示する
- 時定数 $\tau$ を変えたときの減衰曲線を比較する
- ノイズ入りのダミーデータを作る
-
scipy.optimize.curve_fitで指数関数的減衰をフィッティングする - 残差を確認する
- 片対数プロットで指数関数的減衰らしさを見る
という流れで進みます。
図は fig_exp_decay フォルダにPNGとして保存されます。
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.stats import t as student_t
# -----------------------------
# Settings
# -----------------------------
FIG_DIR = "fig_exp_decay"
os.makedirs(FIG_DIR, exist_ok=True)
rng = np.random.default_rng(42)
# Colabで日本語フォントを設定しなくてもよいように、図中の文字は英語にしています。
plt.rcParams["figure.dpi"] = 120
def save_and_show(filename):
path = os.path.join(FIG_DIR, filename)
plt.tight_layout()
plt.savefig(path, dpi=160, bbox_inches="tight")
plt.show()
print(f"saved: {path}")
def exp_decay(t, A, tau, C):
"""
Exponential decay model:
y(t) = C + A * exp(-t / tau)
"""
return C + A * np.exp(-t / tau)
# -----------------------------
# Figure 1: time constant and half-life
# -----------------------------
tau_demo = 2.5
t_demo = np.linspace(0, 8, 400)
y_demo = np.exp(-t_demo / tau_demo)
t_half_demo = tau_demo * np.log(2)
plt.figure(figsize=(8, 5))
plt.plot(t_demo, y_demo, label=r"$e^{-t/\tau}$")
plt.axvline(tau_demo, linestyle="--", label=r"$t=\tau$")
plt.axhline(1 / np.e, linestyle="--", label=r"$1/e \approx 0.368$")
plt.axvline(t_half_demo, linestyle=":", label=r"$t_{1/2}=\tau\ln2$")
plt.axhline(0.5, linestyle=":", label=r"$1/2$")
plt.scatter([tau_demo, t_half_demo], [1 / np.e, 0.5], s=60)
plt.title("Time constant and half-life in exponential decay")
plt.xlabel("time t")
plt.ylabel("relative signal")
plt.ylim(-0.02, 1.05)
plt.legend()
save_and_show("fig1_tau_half_life.png")
# -----------------------------
# Figure 2: different tau values
# -----------------------------
plt.figure(figsize=(8, 5))
for tau in [1.0, 2.5, 5.0]:
plt.plot(t_demo, np.exp(-t_demo / tau), label=fr"$\tau={tau}$")
plt.title("Larger tau means slower decay")
plt.xlabel("time t")
plt.ylabel("relative signal")
plt.ylim(-0.02, 1.05)
plt.legend()
save_and_show("fig2_tau_comparison.png")
# -----------------------------
# Create dummy data
# -----------------------------
A_true = 8.0
tau_true = 2.5
C_true = 1.0
noise_sd = 0.35
t_data = np.linspace(0, 10, 40)
y_clean = exp_decay(t_data, A_true, tau_true, C_true)
y_obs = y_clean + rng.normal(0, noise_sd, size=t_data.size)
# -----------------------------
# Fit the model to the dummy data
# -----------------------------
p0 = [
max(y_obs[0] - y_obs[-1], 0.1), # initial guess for A
2.0, # initial guess for tau
y_obs[-1], # initial guess for C
]
bounds = (
[0, 1e-6, -np.inf], # lower bounds for A, tau, C
[np.inf, np.inf, np.inf], # upper bounds for A, tau, C
)
popt, pcov = curve_fit(
exp_decay,
t_data,
y_obs,
p0=p0,
bounds=bounds,
)
A_hat, tau_hat, C_hat = popt
parameter_se = np.sqrt(np.diag(pcov))
# curve_fitが返す共分散行列から、近似的な95%信頼区間を計算します。
# これは「モデルが正しい」「誤差の仮定がそれなりに妥当」という前提つきの目安です。
dof = len(t_data) - len(popt)
t_value = student_t.ppf(0.975, dof)
ci_A = (A_hat - t_value * parameter_se[0], A_hat + t_value * parameter_se[0])
ci_tau = (tau_hat - t_value * parameter_se[1], tau_hat + t_value * parameter_se[1])
ci_C = (C_hat - t_value * parameter_se[2], C_hat + t_value * parameter_se[2])
t_half_true = tau_true * np.log(2)
t_half_hat = tau_hat * np.log(2)
print("===== True parameters used to create dummy data =====")
print(f"A_true = {A_true:.3f}")
print(f"tau_true = {tau_true:.3f}")
print(f"C_true = {C_true:.3f}")
print(f"t_half_true = {t_half_true:.3f}")
print()
print("===== Estimated parameters from noisy dummy data =====")
print(f"A_hat = {A_hat:.3f} 95% CI [{ci_A[0]:.3f}, {ci_A[1]:.3f}]")
print(f"tau_hat = {tau_hat:.3f} 95% CI [{ci_tau[0]:.3f}, {ci_tau[1]:.3f}]")
print(f"C_hat = {C_hat:.3f} 95% CI [{ci_C[0]:.3f}, {ci_C[1]:.3f}]")
print(f"t_half_hat = {t_half_hat:.3f}")
# -----------------------------
# Figure 3: dummy data and fitted curve
# -----------------------------
t_fit = np.linspace(t_data.min(), t_data.max(), 400)
y_true_fit = exp_decay(t_fit, A_true, tau_true, C_true)
y_fit = exp_decay(t_fit, A_hat, tau_hat, C_hat)
plt.figure(figsize=(8, 5))
plt.scatter(t_data, y_obs, s=35, label="observed dummy data")
plt.plot(t_fit, y_true_fit, linestyle="--", label="true curve")
plt.plot(t_fit, y_fit, label="fitted curve")
plt.axvline(
t_half_hat,
linestyle=":",
label=fr"estimated half-life $\approx {t_half_hat:.2f}$",
)
plt.title("Fitting exponential decay to noisy dummy data")
plt.xlabel("time t")
plt.ylabel("signal y")
plt.legend()
save_and_show("fig3_fit_dummy_data.png")
# -----------------------------
# Figure 4: residuals
# -----------------------------
residuals = y_obs - exp_decay(t_data, A_hat, tau_hat, C_hat)
plt.figure(figsize=(8, 4))
plt.axhline(0, linestyle="--")
plt.scatter(t_data, residuals, s=35)
plt.title("Residuals: observed - fitted")
plt.xlabel("time t")
plt.ylabel("residual")
save_and_show("fig4_residuals.png")
# -----------------------------
# Figure 5: semi-log view
# -----------------------------
# Exponential decay becomes linear after subtracting the baseline and taking log.
# But late noisy points can become unstable when y - C is close to zero.
positive_mask = (y_obs - C_hat) > 0
plt.figure(figsize=(8, 5))
plt.scatter(
t_data[positive_mask],
np.log(y_obs[positive_mask] - C_hat),
s=35,
label=r"$\log(y-C_{\mathrm{fit}})$",
)
plt.plot(
t_fit,
np.log(np.maximum(y_fit - C_hat, 1e-12)),
label="log of fitted exponential part",
)
plt.title("After subtracting baseline, exponential decay becomes linear on a log scale")
plt.xlabel("time t")
plt.ylabel(r"$\log(y-C)$")
plt.legend()
save_and_show("fig5_semilog.png")
実行結果の例
乱数シードを固定しているので、同じ環境ではおおむね次のような結果になるはずです。
===== True parameters used to create dummy data =====
A_true = 8.000
tau_true = 2.500
C_true = 1.000
t_half_true = 1.733
===== Estimated parameters from noisy dummy data =====
A_hat = 7.845 95% CI [7.488, 8.203]
tau_hat = 2.494 95% CI [2.203, 2.785]
C_hat = 1.057 95% CI [0.829, 1.285]
t_half_hat = 1.728
今回は、ダミーデータを作るときに真の値として、
$$
A=8.0,\quad \tau=2.5,\quad C=1.0
$$
を使いました。
フィッティング結果では、
$$
\hat{\tau} \approx 2.494
$$
となっており、真の値 $\tau=2.5$ にかなり近い値が推定できています。
また、半減期は、
$$
t_{1/2} = \tau \log 2
$$
なので、推定された半減期は、
$$
\hat{t}_{1/2} \approx 1.728
$$
です。
真の半減期は、
$$
2.5 \times \log 2 \approx 1.733
$$
なので、こちらもかなり近い値になっています。
図1：時定数と半減期の関係
図1では、シンプルな指数関数的減衰
$$
e^{-t/\tau}
$$
を描いています。
点線で示したように、$t=\tau$ のとき、値は $1/e$ になります。
$$
\frac{1}{e} \approx 0.368
$$
つまり、時定数 $\tau$ は「半分になる時間」ではなく、「約36.8%まで減る時間」です。
一方、半減期 $t_{1/2}$ は、値が $1/2$ になる時間です。
図を見ると、半減期のほうが時定数より少し左にあります。これは、
$$
t_{1/2} = \tau \log 2 \approx 0.693\tau
$$
だからです。
図2：時定数が大きいほど、減衰は遅くなる
図2では、$\tau=1.0$, $\tau=2.5$, $\tau=5.0$ の3つを比較しています。
$\tau$ が小さいと、曲線は急速に下がります。
$\tau$ が大きいと、曲線はゆっくり下がります。
つまり、時定数 $\tau$ は、
どれくらいゆっくり元に戻るか
を表すパラメータです。
指数関数的減衰を見たときに、まず確認したいのはこの $\tau$ です。
図3：ノイズ入りダミーデータをフィッティングする
図3では、青い点がノイズ入りのダミーデータです。
破線は、ダミーデータを作るときに使った「真の曲線」です。
実データでは真の曲線は分かりません。しかし今回はダミーデータなので、答え合わせのために真の曲線も描いています。
実線は、curve_fit によって推定されたフィッティング曲線です。
青い点にはノイズが入っていますが、全体としては指数関数的減衰の形をよく捉えています。
今回のモデルは、
$$
y(t) = C + A e^{-t/\tau}
$$
でした。
フィッティングでは、観測データに合うように $A$, $\tau$, $C$ を推定しています。
フィッティングとは何をしているのか
フィッティングとは、ざっくり言うと、
データに合うように、モデルのパラメータを調整すること
です。
今回の場合、モデルは次の式です。
$$
y(t) = C + A e^{-t/\tau}
$$
この式の中で、未知なのは $A$, $\tau$, $C$ です。
そこでPythonに、
この式が観測データにできるだけ合うように、$A$, $\tau$, $C$ を探してください
と頼んでいます。
curve_fit は、基本的には「観測値」と「モデル予測値」のズレが小さくなるようにパラメータを探します。
このズレを残差と呼びます。
$$
\text{residual}_i = y_i - \hat{y}_i
$$
ここで、$y_i$ は観測値、$\hat{y}_i$ はモデルによる予測値です。
図4：残差を見る
図4は、残差を時間に対してプロットしたものです。
残差は、
$$
\text{観測値} - \text{フィッティング値}
$$
です。
よいフィッティングでは、残差がゼロのまわりにランダムに散らばります。
今回の図では、残差はおおむねゼロのまわりに散らばっており、大きな系統的パターンは見えません。
一方で、もし残差が次のような形をしていたら注意が必要です。
- 前半はずっと正で、後半はずっと負
- 波のような周期的パターンがある
- 時間が経つほど残差が大きくなる
- 一部に極端な外れ値がある
そのような場合、モデルが合っていない、ノイズの性質が変わっている、外れ値がある、などの可能性があります。
図5：指数関数的減衰は、片対数で直線になる
指数関数的減衰は、ベースライン $C$ を引いたあとに対数を取ると、直線になります。
まず、
$$
y(t) - C = A e^{-t/\tau}
$$
です。
両辺の自然対数を取ると、
$$
\log(y(t)-C) = \log A - \frac{t}{\tau}
$$
となります。
これは、
$$
\text{縦軸} = \log(y-C),\quad \text{横軸} = t
$$
とすると、傾きが $-1/\tau$ の直線になることを意味します。
したがって、指数関数的減衰かどうかをざっくり確認するには、片対数プロットが便利です。
ただし、注意点もあります。
時間が十分に経つと、$y-C$ はとても小さくなります。すると、少しのノイズやベースライン推定のズレが、$\log(y-C)$ では大きく見えてしまいます。
図5でも、後半の点は前半より散らばりが大きく見えます。
そのため、片対数プロットは便利ですが、常に万能というわけではありません。
なぜ最初から対数を取って直線フィットしないのか
指数関数的減衰は、対数を取ると直線になるので、
では最初から $\log(y)$ を取って直線フィットすればよいのでは？
と思うかもしれません。
状況によっては、それでもよい場合があります。
ただし、実データでは注意が必要です。
まず、ベースライン $C$ がある場合、対数を取るべきなのは $y$ ではなく $y-C$ です。
$$
\log(y-C)
$$
しかし、$C$ が未知の場合には、先に $C$ を推定する必要があります。
さらに、ノイズがあると、$y-C$ がゼロ以下になる点が出ることもあります。その場合、対数を取ることができません。
また、対数変換をすると、誤差の見え方も変わります。
そのため、今回は
$$
y(t) = C + A e^{-t/\tau}
$$
という元の形のまま、非線形フィッティングを行いました。
実データで注意したいこと
この記事では、仕組みを理解するために、きれいなダミーデータを使いました。
実データでは、もう少し注意が必要です。
1. 本当に単一指数関数でよいとは限らない
現実のデータでは、1つの指数関数だけでは説明できないことがあります。
たとえば、速い成分と遅い成分が混ざっている場合は、
$$
y(t) = C + A_1 e^{-t/\tau_1} + A_2 e^{-t/\tau_2}
$$
のような形が必要になることもあります。
これを二相性減衰、または二重指数関数的減衰として扱う場合があります。
ただし、パラメータが増えるほど、フィッティングは不安定になりやすくなります。
2. ベースラインの推定は重要
指数関数的減衰では、ベースライン $C$ の推定がかなり重要です。
$C$ が少しズレるだけで、$\tau$ や半減期の推定にも影響することがあります。
特に、後半のデータが少ない場合や、十分にベースラインへ近づいていない場合には注意が必要です。
3. フィッティング範囲によって結果が変わる
どの時間範囲を使ってフィッティングするかによって、推定される $\tau$ が変わることがあります。
たとえば、初期にアーチファクトがある場合や、後半に別の遅い成分が混ざる場合、全範囲を単純にフィットすると解釈が難しくなります。
4. 外れ値に注意する
外れ値があると、フィッティング結果が引っ張られることがあります。
残差プロットを見ると、外れ値や系統的なズレに気づきやすくなります。
5. 信頼区間は前提つきの目安
今回のコードでは、curve_fit が返す共分散行列から、近似的な95%信頼区間を計算しました。
ただし、これは、
- モデルが妥当である
- ノイズの仮定が大きく外れていない
- パラメータ推定が局所的に安定している
といった前提に依存します。
実データでは、Bootstrapなどで推定の安定性を調べるのも有用です。
まとめ
この記事では、指数関数的減衰をPythonで確認しました。
重要なポイントは次の通りです。
- 指数関数的減衰は、残っている量に比例して減っていく現象である
- 基本モデルは $y(t)=C + A e^{-t/\tau}$ と書ける
- 時定数 $\tau$ は、ベースラインからの差が約36.8%まで減る時間である
- 半減期 $t_{1/2}$ は、ベースラインからの差が半分になる時間である
- 半減期と時定数の関係は $t_{1/2}=\tau\log2$ である
- ノイズ入りデータでも、適切なモデルを使えば $\tau$ や半減期を推定できる
- フィッティング後は、曲線だけでなく残差も確認するとよい
- 片対数プロットは便利だが、ベースラインやノイズに注意が必要である
指数関数的減衰は、非常に基本的ですが、実験データ解析や数理モデルの理解において重要な形です。
単に式を覚えるだけでなく、図を描き、ダミーデータを作り、フィッティングしてみると、意味がかなり掴みやすくなります。
おまけ：試してみると理解が深まる変更
以下を変更して、図やフィッティング結果がどう変わるか試してみてください。
1. ノイズを大きくする
noise_sd = 0.8
ノイズが大きくなると、$\tau$ の推定や信頼区間がどう変わるでしょうか。
2. 時定数を大きくする
tau_true = 5.0
減衰が遅くなったとき、同じ観測時間 $0 \leq t \leq 10$ で十分に推定できるでしょうか。
3. 観測時間を短くする
t_data = np.linspace(0, 4, 40)
ベースラインに十分近づかないデータでは、$C$ や $\tau$ の推定はどうなるでしょうか。
4. ベースラインをゼロにする
C_true = 0.0
ベースラインがゼロの場合、片対数プロットは見やすくなるでしょうか。
5. 外れ値を入れてみる
たとえば、次のように一部の点だけ大きくずらしてみます。
y_obs[10] += 3.0
残差プロットにどのような変化が出るでしょうか。