0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@ynakayama(Yasuhiro Nakayama)

Euler-Maruyama 法で生成した時系列の積分時間スケールと相関を NumPy で検証する

0
Posted at

確率微分方程式から求めた理論値を、生成した時系列からの推定値と比較します。対象は二領域の線形確率過程で、指標は積分時間スケールと領域間相関の 2 つです。麻酔研究の説明に利用された最小モデルを計算例として用います[1]。

結果を解釈する際は、無限時間の理論値、有限区間で積分した理論値、離散時間の理論値、標本からの推定値を分ける必要があります。同じ「時間スケール」でも計算条件が異なれば数値が変わるためです。本稿では、この差を Python と NumPy で再現します。

Python と NumPy の実行環境

検証用の環境を次のように固定します。Linux または macOS のシェルから実行する例です。

項目 設定
Python 3.12
NumPy 2.5.4
外部ライブラリ NumPy のみ
入力 固定シードから生成した仮想時系列
外部認証 不要 
python3.12 -m venv .venv
.venv/bin/python -m pip install numpy==2.5.4
.venv/bin/python experiment.py

以下のプログラムを experiment.py として保存します。乱数生成には numpy.random.default_rng を利用します[2]。領域間相関には numpy.corrcoef を使い[3]、自己相関の計算には numpy.fft.rfft と numpy.fft.irfft を用います[4]。

二領域の確率微分方程式を定義する

二つの領域の活動を $x_1,x_2$ とし、次の確率微分方程式を使用します。

\begin{aligned}
 dx_1 &=(-\alpha x_1+\beta x_2)\,dt+\sqrt{2D}\,dW_1,\\
 dx_2 &=(\beta x_1-\alpha x_2)\,dt+\sqrt{2D}\,dW_2.
\end{aligned}

$\alpha$ は各領域の減衰率、$\beta$ は領域間相互作用、$D$ は雑音強度です。$W_1,W_2$ は独立したウィーナー過程です。定常分布が存在するように $\alpha>\beta\geq0$、$D>0$ とします。

係数行列の固有値は $-\alpha+\beta$ と $-\alpha-\beta$ であり、どちらも負になります。この条件がシミュレーションの出発点です。

Euler-Maruyama 法で時系列を生成する

Euler-Maruyama 法は、確率微分方程式を時間刻み $\Delta t$ ごとに更新する離散化手法です。独立した標準正規乱数 $\xi_{1,n},\xi_{2,n}$ を使うと、更新式は次のとおりです。

\begin{aligned}
x_{1,n+1}&=x_{1,n}+(-\alpha x_{1,n}+\beta x_{2,n})\Delta t+\sqrt{2D\Delta t}\,\xi_{1,n},\\
x_{2,n+1}&=x_{2,n}+(\beta x_{1,n}-\alpha x_{2,n})\Delta t+\sqrt{2D\Delta t}\,\xi_{2,n}.
\end{aligned}

この離散過程が安定する条件は、二つの固有モードの更新係数の絶対値が 1 未満となる $0<\Delta t<2/(\alpha+\beta)$ です。

初期値を両領域とも 0 とし、最初の 50 時間単位を初期過渡区間として捨てます。ここでの時間単位はモデル上の単位であり、秒やミリ秒には対応づけていません。

自己相関から積分時間スケールを推定する

自己相関は、同じ系列の現在の値と過去の値の統計的な関連をラグごとに表したものです。このモデルの定常正規化自己相関は、ラグ $s\geq0$ に対して次の式になります。

C(s)=\frac{(\alpha+\beta)e^{-(\alpha-\beta)s}+(\alpha-\beta)e^{-(\alpha+\beta)s}}{2\alpha}.

無限区間の積分時間スケール $T_\infty$ と、上限 $L$ までの時間スケール $T_L$ を区別します。

\begin{aligned}
T_\infty&=\int_0^\infty C(s)\,ds
=\frac{\alpha^2+\beta^2}{\alpha(\alpha^2-\beta^2)},\\
T_L&=\int_0^L C(s)\,ds\\
&=\frac{\alpha+\beta}{2\alpha}\frac{1-e^{-(\alpha-\beta)L}}{\alpha-\beta}
+\frac{\alpha-\beta}{2\alpha}\frac{1-e^{-(\alpha+\beta)L}}{\alpha+\beta}.
\end{aligned}

標本自己相関は、各ラグ $k$ の共分散を $N-k$ で割り、ラグ 0 の値で正規化します。平均には標本平均を使います。ゼロ埋めした高速フーリエ変換で自己共分散を求め、上限 $L=5$ まで台形則で積分します。

\hat T_L=\Delta t\left[\frac12+\sum_{k=1}^{K-1}\hat C(k\Delta t)+\frac12\hat C(K\Delta t)\right],\qquad K=L/\Delta t.

標本から得る $\hat T_L$ は、無限区間の $T_\infty$ とも有限区間の $T_L$ とも異なる量です。平均との差し引きや有限標本の変動が含まれます。

定常状態の領域間相関を計算する

二領域の定常共分散行列を $\Sigma$ とすると、連続時間モデルの理論相関は次のとおりです。

\Sigma=\frac{D}{\alpha^2-\beta^2}
\begin{pmatrix}\alpha&\beta\\\beta&\alpha\end{pmatrix},
\qquad S=\frac{\Sigma_{12}}{\sqrt{\Sigma_{11}\Sigma_{22}}}=\frac{\beta}{\alpha}.

時系列からは numpy.corrcoef(sample, rowvar=False)[0, 1] により標本相関 $\hat S$ を求めます。

離散時間の Euler-Maruyama 法にも定常相関が存在します。固有モードの更新係数と定常分散を

\begin{aligned}
\rho_\pm&=1-(\alpha\mp\beta)\Delta t,\\
v_\pm&=\frac{2D\Delta t}{1-\rho_\pm^2}
\end{aligned}

とすると、離散過程の理論相関 $S_E$ とラグ $k$ の自己相関 $C_E(k)$ は

S_E=\frac{v_+-v_-}{v_++v_-},\qquad
C_E(k)=\frac{v_+\rho_+^k+v_-\rho_-^k}{v_++v_-}

です。$C_E(k)$ を同じ台形則で積分した値を $T_E$ とします。これにより、連続時間の理論値との差と、有限標本からの推定誤差を分けて評価できます。

覚醒と麻酔の仮想条件で理論値と推定値を比較する

モデルのパラメーターは原典と同じ 2 条件を使用し、$D=1$ とします。各条件で $\Delta t=0.02,0.01$、観測期間を 200、500 とし、各組み合わせを 20 回試します。乱数シードは 20261024 を基準に、条件と反復番号ごとに変更します。

条件 $\alpha$ $\beta$ $T_\infty$ $T_L$（$L=5$） $S$
覚醒に対応する仮想条件 2.0 1.0 0.8333 0.8283 0.5000
麻酔に対応する仮想条件 4.0 0.4 0.2551 0.2551 0.1000

実験コード

import sys

import numpy as np

SEED = 20261024
D = 1.0
BURN_IN = 50.0
L = 5.0
REPEATS = 20
CONDITIONS = (("覚醒", 2.0, 1.0), ("麻酔", 4.0, 0.4))
STEP_SIZES = (0.02, 0.01)
DURATIONS = (200.0, 500.0)


def continuous_reference(alpha, beta, limit):
    slow = alpha - beta
    fast = alpha + beta
    t_inf = (alpha * alpha + beta * beta) / (alpha * (alpha * alpha - beta * beta))
    t_finite = (
        (fast / (2.0 * alpha)) * (-np.expm1(-slow * limit)) / slow
        + (slow / (2.0 * alpha)) * (-np.expm1(-fast * limit)) / fast
    )
    return t_inf, t_finite, beta / alpha


def euler_reference(alpha, beta, dt, limit, diffusion):
    rho_slow = 1.0 - (alpha - beta) * dt
    rho_fast = 1.0 - (alpha + beta) * dt
    assert abs(rho_slow) < 1.0 and abs(rho_fast) < 1.0
    var_slow = 2.0 * diffusion * dt / (1.0 - rho_slow ** 2)
    var_fast = 2.0 * diffusion * dt / (1.0 - rho_fast ** 2)
    lag = np.arange(round(limit / dt) + 1)
    corr = (var_slow * rho_slow ** lag + var_fast * rho_fast ** lag) / (var_slow + var_fast)
    t_ref = np.trapezoid(corr, dx=dt)
    s_ref = (var_slow - var_fast) / (var_slow + var_fast)
    return t_ref, s_ref


def simulate(alpha, beta, dt, duration, seed, diffusion):
    burn = round(BURN_IN / dt)
    count = round(duration / dt)
    rng = np.random.default_rng(seed)
    noise = rng.standard_normal((burn + count, 2))
    result = np.empty((count, 2))
    x1 = x2 = 0.0
    sigma = np.sqrt(2.0 * diffusion * dt)
    for i, (z1, z2) in enumerate(noise):
        next_x1 = x1 + dt * (-alpha * x1 + beta * x2) + sigma * z1
        next_x2 = x2 + dt * (beta * x1 - alpha * x2) + sigma * z2
        x1, x2 = next_x1, next_x2
        if i >= burn:
            result[i - burn, 0] = x1
            result[i - burn, 1] = x2
    return result


def estimate_time_scale(series, dt, limit):
    centered = series - series.mean()
    n = len(centered)
    max_lag = round(limit / dt)
    assert max_lag < n
    freq = np.fft.rfft(centered, n=2*n)
    autocov = np.fft.irfft(freq * freq.conj(), n=2*n)[:max_lag+1].real
    autocov /= n - np.arange(max_lag + 1)
    acf = autocov / autocov[0]
    return float(np.trapezoid(acf, dx=dt))


def main():
    print("Python", ".".join(map(str, sys.version_info[:3])), "NumPy", np.__version__)
    for condition_i, (label, alpha, beta) in enumerate(CONDITIONS):
        t_inf, t_lim, s_cont = continuous_reference(alpha, beta, L)
        print(f"{label}: T∞={t_inf:.4f} TL={t_lim:.4f} S={s_cont:.4f}")
        print("dt   期間   TE     平均T   MAE_T  SE     平均S   MAE_S")
        for dt_i, dt in enumerate(STEP_SIZES):
            for duration_i, duration in enumerate(DURATIONS):
                t_euler, s_euler = euler_reference(alpha, beta, dt, L, D)
                times, cors = [], []
                for repeat in range(REPEATS):
                    seed = SEED + condition_i * 10000 + dt_i * 1000 + duration_i * 100 + repeat
                    sample = simulate(alpha, beta, dt, duration, seed, D)
                    times.append(estimate_time_scale(sample[:, 0], dt, L))
                    cors.append(float(np.corrcoef(sample, rowvar=False)[0, 1]))
                t = np.array(times)
                s = np.array(cors)
                print(f"{dt:.2f} {duration:.0f} {t_euler:.4f} {t.mean():.4f} {np.abs(t - t_euler).mean():.4f} {s_euler:.4f} {s.mean():.4f} {np.abs(s - s_euler).mean():.4f}")


if __name__ == "__main__":
    main()

TE と SE は離散時間の理論値、平均T と 平均S は 20 回の推定値の算術平均です。MAE_T と MAE_S は、各反復と離散時間の理論値の絶対差を 20 回分平均した値です。

実行例

次の数値は Python 3.13.5、NumPy 2.3.5 で実際に実行した結果です。冒頭に示した Python 3.12、NumPy 2.5.4 の組み合わせでは未実行であり、その環境の結果としては提示しません。表示値は小数第 4 位に丸めています。

Python 3.13.5 NumPy 2.3.5
覚醒: T∞=0.8333 TL=0.8283 S=0.5000
dt   期間   TE     平均T   MAE_T  SE     平均S   MAE_S
0.02 200 0.8160 0.7968 0.1648 0.4923 0.5000 0.0379
0.02 500 0.8160 0.8225 0.1015 0.4923 0.4958 0.0195
0.01 200 0.8222 0.8059 0.0977 0.4962 0.4852 0.0294
0.01 500 0.8222 0.7474 0.1212 0.4962 0.4893 0.0191
麻酔: T∞=0.2551 TL=0.2551 S=0.1000
dt   期間   TE     平均T   MAE_T  SE     平均S   MAE_S
0.02 200 0.2449 0.2426 0.0617 0.0959 0.1000 0.0281
0.02 500 0.2449 0.2575 0.0465 0.0959 0.0934 0.0174
0.01 200 0.2500 0.2104 0.0532 0.0980 0.0869 0.0251
0.01 500 0.2500 0.2240 0.0441 0.0980 0.0987 0.0227

観測期間と時間刻みを変更して推定誤差を測る

時間刻みを 0.02 から 0.01 にすると、両条件で $T_E$ は有限区間の理論値 $T_L$ に、$S_E$ は連続時間の理論相関 $S$ に近づいています。例えば覚醒の仮想条件では $T_E$ が 0.8160 から 0.8222 に、$S_E$ が 0.4923 から 0.4962 へ変わります。これは数値更新と台形積分を組み合わせた離散計算の基準値の差です。

観測期間を 200 から 500 に増やした場合、領域間相関の平均絶対誤差は各条件・時間刻みで小さくなっています。一方、時間スケールの平均絶対誤差は覚醒条件の $\Delta t=0.01$ で 0.0977 から 0.1212 に増えています。20 回の実験で誤差が必ず単調減少するとは限らず、自己相関の長いラグは標本変動の影響を受けます。

固定した乱数シードは再実行による比較を可能にしますが、別の乱数シードや反復数では推定値の集計結果が変わります。観測期間による誤差の傾向をより厳密に調べる場合は、反復数を増やし、平均絶対誤差の変動も確認する必要があります。

有限観測と離散化による誤差を区別する

時間スケールについて、比較対象を次のように整理できます。

比較 意味
$T_\infty-T_L$ 積分範囲を $L$ で打ち切った差
$T_L-T_E$ Euler 離散化と離散時刻の台形積分による基準値の差
$\hat T_L-T_E$ 有限長時系列からの推定値と離散過程の理論値の差
$S-S_E$ 連続時間と Euler 離散過程の相関の差
$\hat S-S_E$ 標本相関の推定誤差

例えば覚醒条件の $\Delta t=0.02$ では、$T_\infty=0.8333$、$T_L=0.8283$、$T_E=0.8160$ です。最初の差約 0.0050 は積分範囲の打ち切り、次の差約 0.0123 は離散更新と求積法の差です。標本推定値のばらつきはさらに別に評価します。

また、標本自己相関は平均を標本から推定し、各ラグで $N-k$ によって補正しています。この推定方法がもたらす有限標本の偏りも $\hat T_L-T_E$ に含まれます。MAE_T はこの差の絶対値を反復平均したものであり、単なる離散化誤差ではありません。

まとめ

今回の二領域モデルでは、積分時間スケールと領域間相関の連続時間の理論値を解析的に求められます。しかし、実装で生成するのは時間刻みを持つ離散過程であり、観測できるのは有限長の時系列です。

比較を正しく行うには、まず $T_\infty$ と $T_L$ を区別し、次に $T_E$ と $S_E$ を離散過程の基準値として求めます。そのうえで $\hat T_L$ と $\hat S$ の反復間のばらつきを評価することで、積分範囲の打ち切り、離散化、有限標本の影響を混同せずに検証できます。

この仮想モデルの数値は神経活動の実測値を再現するものではありません。時系列の推定方法を評価するとき、生成モデルの理論値と観測手続きを対応づけるための計算例です。

参考文献

  1. id774, 麻酔で共通する神経活動だけでは感じ方を特定できない（2026-10-15）. https://blog.id774.net/entry/2026/10/15/5764/
  2. NumPy Developers, Random Generator. https://numpy.org/doc/stable/reference/random/generator.html
  3. NumPy Developers, numpy.corrcoef. https://numpy.org/doc/stable/reference/generated/numpy.corrcoef.html
  4. NumPy Developers, Discrete Fourier Transform. https://numpy.org/doc/stable/reference/routines.fft.html
0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?