不等間隔の測定値を台形則で積分するには、trapezoid(y, x=t)のように実際の時刻を渡す。xを省略すると間隔1とみなされる。また、この関数は時刻を自動で並べ替えない。時刻と値の対応、順序、単位、欠測区間を確認してから積分しよう。
最小例で確かめる
以下は説明用に作った小さなデータである。実測データや実行速度の測定結果ではない。コード全体をexample.pyとして保存すれば、入力ファイルを別途用意せずに実行できる。assertは、この例で成り立つべき形や値を確認するために入れてある。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid
def integrate_observed(t, y):
t, y = np.asarray(t, dtype=float), np.asarray(y, dtype=float)
if t.ndim != 1 or y.shape != t.shape or t.size < 2:
raise ValueError("matching one-dimensional arrays are required")
if not (np.isfinite(t).all() and np.isfinite(y).all()):
raise ValueError("missing or infinite sample")
if not (np.diff(t) > 0).all():
raise ValueError("time must increase strictly")
return trapezoid(y, x=t)
t = np.array([0., 0.5, 2., 3.]) # s
y = 2 * t + 1 # 任意の流量単位
area = integrate_observed(t, y)
print("with timestamps:", float(area))
print("unit spacing (wrong here):", float(trapezoid(y)))
print("reversed path:", float(trapezoid(y[::-1], x=t[::-1])))
for label, tx, yy in [("unordered", t[[0,2,1,3]], y[[0,2,1,3]]),
("missing", t, [1., 2., np.nan, 7.])]:
try:
integrate_observed(tx, yy)
except ValueError:
print("expected rejection:", label)
else:
raise AssertionError("invalid sample was accepted")
np.testing.assert_allclose(area, 12., rtol=0, atol=1e-12)
fig, ax = plt.subplots(figsize=(7, 4), layout="constrained")
ax.plot(t, y, "o-", color="#a83262", label="synthetic samples")
ax.fill_between(t, y, alpha=0.18, color="#a83262")
ax.vlines(t, 0, y, color="0.5", linewidth=1)
ax.set(xlabel="Time [s]", ylabel="Flow [arb. unit / s]", ylim=(0, 8),
title="Unequal time intervals: trapezoid area = 12")
ax.legend()
fig.savefig("PY082-trapezoid.png", dpi=160)
plt.close(fig)
実行結果
with timestamps: 12.0
unit spacing (wrong here): 11.0
reversed path: -12.0
expected rejection: unordered
expected rejection: missing
幅が異なる台形を足す
台形則は、隣り合う測定値の平均にその区間の幅を掛け、各区間について足し合わせる。ここでは幅が0.5、1.5、1.0秒なので、どの区間も同じ幅だとする計算は合わない。yが流量、tが秒なら、積分結果の単位は流量×秒になる。
例のy=2t+1は説明用に作った直線であり、区分的に直線で結ぶ台形則と一致するため、積分値12を解析的に確認できる。一般の曲線や実測データで同じ精度を保証する例ではない。表示した図の面積は、測定点を直線で結ぶという仮定に基づくものだと理解してほしい。
trapezoidは入力順に進む
xを渡した場合、trapezoidは与えた順番の差分を使う。時刻と値を両方逆順にすると、経路を逆向きに積分したことになり、結果は-12になる。一般の曲線積分ではこの挙動が必要だが、通常の時系列の積算量を求めたいなら、時間が厳密に増加しているか確認する。
ログが単に順不同なら、時刻と測定値を同じ並べ替えindexで並べる。一方、装置の時計がリセットされた、複数の走行記録が混ざった、時差が混在する場合には、ソートするだけでは本来つながらないデータを接続してしまう。時刻が戻った理由を調べ、別の測定区間に分ける必要がある。
欠測を削ると空白区間を橋渡しする
欠測点を落として残った両端を積分すると、その間を直線で結んだことになる。長い欠測区間の中で値がどう変化したかは観測されていないので、その面積を実測から分かった値として扱うべきではない。例の関数は非有限値を拒否し、判断を後段へ持ち越さない設計にした。
欠測を許す処理なら、連続して観測できた区間ごとに積分し、積分した時間範囲も報告する方法がある。一定時間以上の空白を区切る、物理モデルに基づいて補間するなど、採用する規則を明示しよう。何秒の欠測まで許すかは、現象の変化速度や測定目的による。
数値誤差と測定誤差を分ける
標本点の間に急な変化があれば、台形則の近似誤差が大きくなる。配列を細かく補間して点を増やすだけでは、新しい観測情報は増えない。積分法を高度にしたことで未知のピークを回復できた、と考えないこと。測定間隔そのものが十分かも検討する必要がある。
時刻の重複を許すかも確認する。幅0の区間は面積に寄与しなくても、同じ時刻に値が二つある理由は未解決である。ここでは一つの時系列として厳密増加を要求した。多チャンネルの積分へ広げる場合には、時刻方向のaxisを明示し、出力shapeと各列の単位を合わせて確認するとよい。
動作確認環境と参考資料
Linux・CPython 3.12.14、NumPy 2.3.5、pandas 2.2.3、SciPy 1.17.0、Matplotlib 3.10.8の環境で掲載コードを実行した。使用するライブラリはコード冒頭のimportを参照してほしい。公式資料の最新版と、この実行確認版は区別している。数値の末尾や表の表示幅は環境によって変わることがある。