確率的に変動する量を時間の経過とともに追跡したい場面は、工学・物理学・金融・生物学など幅広い分野で現れます。株価の変動、通信チャネルのノイズ、粒子の拡散運動――これらはすべて 確率過程（stochastic process） という統一的な数学的枠組みで扱うことができます。
確率過程は確率論の中でも特に豊かな構造を持つ分野であり、マルコフ連鎖、ブラウン運動、ポアソン過程などの個別のトピックを学ぶ前に、その全体像を把握しておくことが非常に重要です。本記事では、確率過程の厳密な定義から始めて、分類体系、定常性、エルゴード性、パワースペクトル密度まで丁寧に解説し、主要な確率過程のサンプルパスをPythonで生成・可視化します。
本記事の内容
- 確率過程の厳密な定義 $\{X(t), t \in T\}$
- 分類法（離散/連続時間 × 離散/連続状態）
- 定常過程（強定常/弱定常）の定義と違い
- エルゴード性の定義と意味
- 自己相関関数とパワースペクトル密度（ウィーナー=ヒンチンの定理）
- 主要な確率過程の関係図
- Pythonでの各種確率過程のサンプルパス生成・可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
確率過程とは
直感的な理解
確率過程とは、大雑把に言えば「時間とともにランダムに変化する量」を数学的にモデル化したものです。たとえば、ある株式の価格 $S(t)$ は時刻 $t$ ごとに異なる値をとりますが、明日の株価がいくらになるかは確定的には分かりません。しかし、その統計的な振る舞い（平均リターン、ボラティリティなど）は確率的に記述できます。
イメージとしては、サイコロを毎秒振って出た目を記録するような実験を考えてみてください。1回の実験で得られる時系列データ（サンプルパス）は毎回異なりますが、その背後にある確率法則は同じです。
数学的定義
確率過程を厳密に定義します。
定義（確率過程）: 確率空間 $(\Omega, \mathcal{F}, P)$ とパラメータ集合 $T$ が与えられたとき、各 $t \in T$ に対して定義された確率変数 $X(t) : \Omega \to S$ の族
$$ \{X(t, \omega) : t \in T, \omega \in \Omega\} $$
を 確率過程（stochastic process） と呼びます。ここで、
- $T$ は パラメータ集合（多くの場合、時間を表す）
- $S$ は 状態空間（確率変数がとりうる値の集合）
- $\omega \in \Omega$ は 標本点（ランダムさの源）
この定義には2つの見方があります。
見方1: 時刻を固定 — ある時刻 $t_0$ を固定すると、$X(t_0, \cdot)$ は通常の確率変数になります。つまり、確率過程は「確率変数の族（ファミリー）」です。
見方2: 標本点を固定 — ある $\omega_0 \in \Omega$ を固定すると、$X(\cdot, \omega_0)$ は $T$ から $S$ への写像（関数）になります。これを サンプルパス（sample path） または 実現値（realization） と呼びます。
確率過程の分類
確率過程は、パラメータ集合 $T$ と状態空間 $S$ の性質に応じて4つに分類できます。
パラメータ集合（時間）による分類
- 離散時間過程: $T = \{0, 1, 2, \dots\}$ や $T = \mathbb{Z}$。時刻は飛び飛びの値をとる。$X_n$ と添字で表記することが多い。
- 連続時間過程: $T = [0, \infty)$ や $T = \mathbb{R}$。時刻は連続的な値をとる。$X(t)$ と関数表記することが多い。
状態空間による分類
- 離散状態過程: $S$ が可算集合（例: $S = \{0, 1, 2, \dots\}$, $S = \{A, B, C\}$）。連鎖（chain） と呼ばれることが多い。
- 連続状態過程: $S = \mathbb{R}$ や $S = \mathbb{R}^d$。
4分類の表
| 離散状態 | 連続状態 | |
|---|---|---|
| 離散時間 | ベルヌーイ過程、マルコフ連鎖、ランダムウォーク | 離散時間ガウス過程 |
| 連続時間 | ポアソン過程、連続時間マルコフ連鎖 | ブラウン運動、オルンシュタイン=ウーレンベック過程 |
有限次元分布による特徴づけ
確率過程の確率法則は、すべての有限次元分布の族
$$ F_{t_1, t_2, \dots, t_n}(x_1, x_2, \dots, x_n) = P(X(t_1) \leq x_1, X(t_2) \leq x_2, \dots, X(t_n) \leq x_n) $$
によって完全に特徴づけられます（コルモゴロフの拡張定理）。ここで $n \geq 1$ は任意の正整数、$t_1, t_2, \dots, t_n \in T$ は任意の時点です。
この有限次元分布族は以下の 整合性条件（consistency condition） を満たす必要があります。
対称性条件: 任意の置換 $\sigma$ に対して、
$$ F_{t_{\sigma(1)}, \dots, t_{\sigma(n)}}(x_{\sigma(1)}, \dots, x_{\sigma(n)}) = F_{t_1, \dots, t_n}(x_1, \dots, x_n) $$
周辺化条件: $x_n \to \infty$ のとき、
$$ F_{t_1, \dots, t_{n-1}, t_n}(x_1, \dots, x_{n-1}, \infty) = F_{t_1, \dots, t_{n-1}}(x_1, \dots, x_{n-1}) $$
平均関数と自己相関関数
確率過程の基本的な統計量として、平均関数と自己相関関数を定義します。
平均関数
確率過程 $\{X(t)\}$ の 平均関数（mean function） は、
$$ \mu_X(t) = E[X(t)] $$
として定義されます。各時刻 $t$ での期待値を取ったものです。
自己相関関数
自己相関関数（autocorrelation function） は、2時点の関係を記述する量で、
$$ R_{XX}(t_1, t_2) = E[X(t_1) X(t_2)] $$
と定義されます。
自己共分散関数
自己共分散関数（autocovariance function） は、
$$ C_{XX}(t_1, t_2) = \text{Cov}(X(t_1), X(t_2)) = E[(X(t_1) – \mu_X(t_1))(X(t_2) – \mu_X(t_2))] $$
と定義されます。展開すると、
$$ \begin{align} C_{XX}(t_1, t_2) &= E[X(t_1) X(t_2) – \mu_X(t_1) X(t_2) – X(t_1) \mu_X(t_2) + \mu_X(t_1) \mu_X(t_2)] \\ &= E[X(t_1) X(t_2)] – \mu_X(t_1) E[X(t_2)] – E[X(t_1)] \mu_X(t_2) + \mu_X(t_1) \mu_X(t_2) \\ &= R_{XX}(t_1, t_2) – \mu_X(t_1) \mu_X(t_2) \end{align} $$
という関係が成り立ちます。
定常過程
強定常過程
定義（強定常過程）: 確率過程 $\{X(t)\}$ が 強定常（strictly stationary） であるとは、任意の $n \geq 1$、任意の $t_1, \dots, t_n \in T$、任意のシフト量 $\tau$ に対して、
$$ F_{t_1 + \tau, t_2 + \tau, \dots, t_n + \tau}(x_1, x_2, \dots, x_n) = F_{t_1, t_2, \dots, t_n}(x_1, x_2, \dots, x_n) $$
が成り立つことをいいます。
つまり、確率過程の有限次元分布が 時間シフト $\tau$ に対して不変 であるということです。直感的に言えば、「統計的な性質が時間の経過とともに変化しない」過程です。
強定常過程の性質を整理します。$n = 1$ の場合、
$$ F_{t + \tau}(x) = F_t(x) \quad \forall \tau $$
となるため、$X(t)$ の分布は $t$ に依存しません。したがって平均関数は定数になります。
$$ \mu_X(t) = \mu_X = \text{const.} $$
$n = 2$ の場合、
$$ F_{t_1 + \tau, t_2 + \tau}(x_1, x_2) = F_{t_1, t_2}(x_1, x_2) $$
ここで $\tau = -t_1$ とおくと、
$$ F_{0, t_2 – t_1}(x_1, x_2) = F_{t_1, t_2}(x_1, x_2) $$
となるため、2次元分布は時刻差 $t_2 – t_1$ のみに依存します。したがって自己相関関数も時刻差のみの関数になります。
$$ R_{XX}(t_1, t_2) = R_{XX}(t_2 – t_1) $$
弱定常過程（広義定常過程）
定義（弱定常過程）: 確率過程 $\{X(t)\}$ が 弱定常（wide-sense stationary, WSS） であるとは、以下の2条件を満たすことをいいます。
- 平均が定数: $E[X(t)] = \mu_X = \text{const.}$ （$t$ に依存しない）
- 自己相関関数が時刻差のみに依存: $R_{XX}(t_1, t_2) = R_{XX}(t_2 – t_1) = R_{XX}(\tau)$
ここで $\tau = t_2 – t_1$ は ラグ（lag） と呼ばれます。
強定常 ⇒ 弱定常 は上で示した通り成り立ちます。しかし、逆は一般には成り立ちません。弱定常は1次と2次のモーメントのみの条件であり、高次のモーメントや分布全体の不変性は要求していません。
ただし、ガウス過程（有限次元分布がすべて正規分布である過程）の場合は、正規分布が平均と分散（共分散）で完全に特徴づけられるため、弱定常 ⇔ 強定常 が成り立ちます。
弱定常過程の自己相関関数の性質
弱定常過程の自己相関関数 $R_{XX}(\tau)$ は以下の性質を持ちます。
性質1: 偶関数
$$ R_{XX}(-\tau) = E[X(t)X(t + \tau)] = E[X(t – \tau)X(t)] = R_{XX}(\tau) $$
2番目の等号は弱定常性から従います。
性質2: 原点で最大
$$ |R_{XX}(\tau)| \leq R_{XX}(0) $$
これはコーシー=シュワルツの不等式から導出できます。
$$ |R_{XX}(\tau)|^2 = |E[X(t)X(t+\tau)]|^2 \leq E[|X(t)|^2] \cdot E[|X(t+\tau)|^2] = R_{XX}(0)^2 $$
性質3: 非負定値
任意の $n$、任意の定数 $a_1, \dots, a_n$、任意の $t_1, \dots, t_n$ に対して、
$$ \sum_{i=1}^{n} \sum_{j=1}^{n} a_i a_j R_{XX}(t_i – t_j) \geq 0 $$
これは、$Y = \sum_{i=1}^{n} a_i X(t_i)$ とおくと、$E[Y^2] \geq 0$ から直ちに従います。
エルゴード性
エルゴード性の直感
定常過程では、統計量（平均、分散、相関など）はアンサンブル（集合）平均として定義されます。つまり、同じ確率過程の独立な実現を多数集めて平均を取る操作です。しかし、実際の応用では、1本のサンプルパスしか観測できないことが多いです。
エルゴード性（ergodicity） とは、「1本のサンプルパスの時間平均がアンサンブル平均と一致する」という性質です。エルゴード性が成り立てば、1回の長い観測から統計量を推定できます。
平均に関するエルゴード性
弱定常過程 $\{X(t)\}$ の 時間平均 を
$$ \langle X \rangle_T = \frac{1}{2T} \int_{-T}^{T} X(t) \, dt $$
と定義します。平均に関するエルゴード的（ergodic in the mean） であるとは、
$$ \lim_{T \to \infty} \langle X \rangle_T = \mu_X \quad \text{（確率収束の意味で）} $$
が成り立つことをいいます。
この条件が成り立つための十分条件は、
$$ \lim_{T \to \infty} \frac{1}{2T} \int_{-T}^{T} \left(1 – \frac{|\tau|}{2T}\right) C_{XX}(\tau) \, d\tau = 0 $$
です。ここで $C_{XX}(\tau) = R_{XX}(\tau) – \mu_X^2$ は自己共分散関数です。
特に、$C_{XX}(\tau) \to 0$ ($|\tau| \to \infty$) であれば、上の条件は満たされます。直感的には、「十分離れた時刻の確率変数が無相関になる」過程はエルゴード的です。
相関に関するエルゴード性
同様に、自己相関に関するエルゴード性は、時間平均
$$ \hat{R}_{XX}(\tau) = \lim_{T \to \infty} \frac{1}{2T} \int_{-T}^{T} X(t) X(t + \tau) \, dt $$
がアンサンブル自己相関 $R_{XX}(\tau)$ に確率収束することを要求します。
ウィーナー=ヒンチンの定理
パワースペクトル密度
弱定常過程の周波数領域での特徴づけとして、パワースペクトル密度（power spectral density, PSD） を導入します。
定義: 弱定常過程 $\{X(t)\}$ の自己相関関数 $R_{XX}(\tau)$ が絶対可積分であるとき、そのフーリエ変換
$$ S_{XX}(f) = \int_{-\infty}^{\infty} R_{XX}(\tau) e^{-j2\pi f \tau} \, d\tau $$
を パワースペクトル密度 と呼びます。ここで $f$ は周波数、$j = \sqrt{-1}$ は虚数単位です。
ウィーナー=ヒンチンの定理の導出
定理（ウィーナー=ヒンチン, Wiener-Khinchin）: 弱定常過程 $\{X(t)\}$ のパワースペクトル密度 $S_{XX}(f)$ と自己相関関数 $R_{XX}(\tau)$ はフーリエ変換対をなす。
$$ S_{XX}(f) = \int_{-\infty}^{\infty} R_{XX}(\tau) e^{-j2\pi f \tau} \, d\tau $$
$$ R_{XX}(\tau) = \int_{-\infty}^{\infty} S_{XX}(f) e^{j2\pi f \tau} \, df $$
導出: 確率過程 $X(t)$ の切断フーリエ変換を
$$ X_T(f) = \int_{-T}^{T} X(t) e^{-j2\pi f t} \, dt $$
と定義します。$|X_T(f)|^2$ の期待値を計算します。
$$ \begin{align} E[|X_T(f)|^2] &= E\left[\int_{-T}^{T} X(t) e^{-j2\pi f t} \, dt \cdot \int_{-T}^{T} X(s)^* e^{j2\pi f s} \, ds\right] \\ &= \int_{-T}^{T} \int_{-T}^{T} E[X(t) X(s)^*] e^{-j2\pi f(t-s)} \, dt \, ds \\ &= \int_{-T}^{T} \int_{-T}^{T} R_{XX}(t – s) e^{-j2\pi f(t-s)} \, dt \, ds \end{align} $$
最後の等号で弱定常性 $E[X(t)X(s)^*] = R_{XX}(t – s)$ を用いました。$\tau = t – s$ と変数変換すると、
$$ E[|X_T(f)|^2] = \int_{-2T}^{2T} \left(2T – |\tau|\right) R_{XX}(\tau) e^{-j2\pi f \tau} \, d\tau $$
両辺を $2T$ で割って $T \to \infty$ の極限を取ると、
$$ \lim_{T \to \infty} \frac{E[|X_T(f)|^2]}{2T} = \int_{-\infty}^{\infty} R_{XX}(\tau) e^{-j2\pi f \tau} \, d\tau = S_{XX}(f) $$
左辺は単位時間あたりの平均パワーの周波数分布を表しており、これが $R_{XX}(\tau)$ のフーリエ変換に等しいことが示されました。
PSDの性質:
- $S_{XX}(f) \geq 0$（非負性）: パワーは非負
- $S_{XX}(-f) = S_{XX}(f)$（偶関数性）: $R_{XX}(\tau)$ が実偶関数であることから従う
- $\int_{-\infty}^{\infty} S_{XX}(f) \, df = R_{XX}(0) = E[|X(t)|^2]$（全パワー）
主要な確率過程の関係
ここで、主要な確率過程の関係を整理します。
ベルヌーイ過程
独立同分布な $\{0, 1\}$ 値確率変数列 $\{X_n\}$ で、$P(X_n = 1) = p$ です。最も単純な離散時間・離散状態の確率過程です。
ランダムウォーク
ベルヌーイ過程の累積和として構成されます。$S_n = \sum_{k=1}^{n} Y_k$ ここで $Y_k$ は独立同分布で $P(Y_k = +1) = p$, $P(Y_k = -1) = 1 – p$ です。離散時間・離散状態の過程で、マルコフ性を持ちます。
マルコフ過程
未来の状態が現在の状態のみに依存し、過去の履歴には依存しない過程です。数学的には、
$$ P(X(t_{n+1}) \leq x \mid X(t_n) = x_n, \dots, X(t_1) = x_1) = P(X(t_{n+1}) \leq x \mid X(t_n) = x_n) $$
ランダムウォークはマルコフ過程の特殊例です。
ポアソン過程
強度 $\lambda > 0$ のポアソン過程 $\{N(t)\}$ は、$[0, t]$ の間に起こるイベントの回数を数える連続時間・離散状態の過程です。独立増分とポアソン分布 $N(t) \sim \text{Poi}(\lambda t)$ を持ちます。
ガウス過程
任意の有限次元分布が多次元正規分布に従う過程です。平均関数 $\mu(t)$ と共分散関数 $C(t_1, t_2)$ で完全に特徴づけられます。
ブラウン運動（ウィーナー過程）
連続時間・連続状態のガウス過程で、独立定常増分を持ちます。ランダムウォークの連続時間極限として得られます（Donskerの不変原理）。物理学、金融工学、確率解析の基盤です。
関係図
ベルヌーイ過程 ──(累積和)──→ ランダムウォーク ──(連続極限)──→ ブラウン運動
│ │
├─── マルコフ過程 ←──────────┤
│ │
ポアソン過程 ←──(計数過程)──── 独立増分過程 ─────────────────────┘
│
ガウス過程 ←┘
Pythonでの各種確率過程のサンプルパス生成
各種確率過程のサンプルパスを生成し、可視化します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
fig, axes = plt.subplots(3, 2, figsize=(14, 12))
# --- (1) ベルヌーイ過程 ---
n_steps = 100
p = 0.5
n_paths = 5
for i in range(n_paths):
bernoulli_path = np.random.binomial(1, p, size=n_steps)
axes[0, 0].step(range(n_steps), bernoulli_path, alpha=0.7, linewidth=0.8)
axes[0, 0].set_title("Bernoulli Process (p=0.5)")
axes[0, 0].set_xlabel("Time step n")
axes[0, 0].set_ylabel("X_n")
axes[0, 0].set_ylim(-0.2, 1.2)
# --- (2) ランダムウォーク ---
n_steps = 200
for i in range(n_paths):
steps = np.random.choice([-1, 1], size=n_steps)
random_walk = np.cumsum(steps)
random_walk = np.insert(random_walk, 0, 0)
axes[0, 1].plot(range(n_steps + 1), random_walk, alpha=0.7, linewidth=0.8)
axes[0, 1].set_title("Random Walk (symmetric)")
axes[0, 1].set_xlabel("Time step n")
axes[0, 1].set_ylabel("S_n")
axes[0, 1].axhline(y=0, color='k', linestyle='--', linewidth=0.5)
# --- (3) ポアソン過程 ---
T_max = 10.0
lam = 3.0 # 強度パラメータ
for i in range(n_paths):
# 指数分布に従う到着間隔を生成
n_events = np.random.poisson(lam * T_max * 2)
inter_arrivals = np.random.exponential(1.0 / lam, size=n_events)
arrival_times = np.cumsum(inter_arrivals)
arrival_times = arrival_times[arrival_times <= T_max]
# 階段関数として描画
t_plot = np.concatenate([[0], np.repeat(arrival_times, 2), [T_max]])
n_plot = np.repeat(range(len(arrival_times) + 1), 2)
axes[1, 0].plot(t_plot, n_plot, alpha=0.7, linewidth=0.8)
axes[1, 0].set_title(f"Poisson Process (λ={lam})")
axes[1, 0].set_xlabel("Time t")
axes[1, 0].set_ylabel("N(t)")
# --- (4) マルコフ連鎖（3状態） ---
# 遷移確率行列
P_mat = np.array([
[0.7, 0.2, 0.1],
[0.3, 0.4, 0.3],
[0.1, 0.3, 0.6]
])
states = [0, 1, 2]
n_steps_mc = 100
for i in range(n_paths):
path = [np.random.choice(states)]
for _ in range(n_steps_mc - 1):
current = path[-1]
next_state = np.random.choice(states, p=P_mat[current])
path.append(next_state)
axes[1, 1].step(range(n_steps_mc), path, alpha=0.7, linewidth=0.8)
axes[1, 1].set_title("Markov Chain (3 states)")
axes[1, 1].set_xlabel("Time step n")
axes[1, 1].set_ylabel("State")
axes[1, 1].set_yticks([0, 1, 2])
# --- (5) ブラウン運動 ---
T_bm = 1.0
n_points = 1000
dt = T_bm / n_points
t_bm = np.linspace(0, T_bm, n_points + 1)
for i in range(n_paths):
dW = np.random.normal(0, np.sqrt(dt), size=n_points)
W = np.cumsum(dW)
W = np.insert(W, 0, 0)
axes[2, 0].plot(t_bm, W, alpha=0.7, linewidth=0.8)
axes[2, 0].set_title("Brownian Motion (Wiener Process)")
axes[2, 0].set_xlabel("Time t")
axes[2, 0].set_ylabel("W(t)")
axes[2, 0].axhline(y=0, color='k', linestyle='--', linewidth=0.5)
# --- (6) オルンシュタイン=ウーレンベック過程 ---
# dX = theta*(mu - X)dt + sigma*dW
theta = 5.0
mu_ou = 0.0
sigma_ou = 1.0
for i in range(n_paths):
X_ou = np.zeros(n_points + 1)
X_ou[0] = np.random.normal(0, 1) # 初期値
for k in range(n_points):
dW_k = np.random.normal(0, np.sqrt(dt))
X_ou[k + 1] = X_ou[k] + theta * (mu_ou - X_ou[k]) * dt + sigma_ou * dW_k
axes[2, 1].plot(t_bm, X_ou, alpha=0.7, linewidth=0.8)
axes[2, 1].set_title("Ornstein-Uhlenbeck Process")
axes[2, 1].set_xlabel("Time t")
axes[2, 1].set_ylabel("X(t)")
axes[2, 1].axhline(y=mu_ou, color='k', linestyle='--', linewidth=0.5)
plt.tight_layout()
plt.savefig("stochastic_processes_sample_paths.png", dpi=150, bbox_inches="tight")
plt.show()
このコードでは、6つの主要な確率過程のサンプルパスを生成しています。各過程の特徴を観察してみましょう。
- ベルヌーイ過程: 各時刻で0または1の値をとり、時刻間で独立
- ランダムウォーク: 累積和なので連続的な変動を示し、原点から離れていく傾向がある
- ポアソン過程: 階段関数で、ジャンプの間隔が指数分布に従う
- マルコフ連鎖: 離散状態間を遷移し、次の状態は現在の状態のみに依存
- ブラウン運動: 連続だが至る所微分不可能で、不規則な動き
- オルンシュタイン=ウーレンベック過程: 平均回帰性を持ち、$\mu = 0$ の周りを振動
弱定常過程の自己相関関数とPSDの可視化
弱定常過程の例として、自己相関関数とパワースペクトル密度のフーリエ変換対の関係を可視化します。
import numpy as np
import matplotlib.pyplot as plt
# 指数減衰型の自己相関関数を持つ弱定常過程
# R_XX(tau) = sigma^2 * exp(-alpha * |tau|)
# この場合のPSD: S_XX(f) = 2 * sigma^2 * alpha / (alpha^2 + (2*pi*f)^2)
sigma2 = 1.0 # 分散
alpha = 2.0 # 減衰定数
# 自己相関関数
tau = np.linspace(-5, 5, 1000)
R_XX = sigma2 * np.exp(-alpha * np.abs(tau))
# パワースペクトル密度（解析的な式）
f = np.linspace(-5, 5, 1000)
S_XX = 2 * sigma2 * alpha / (alpha**2 + (2 * np.pi * f)**2)
# パワースペクトル密度（FFTによる数値計算）
dtau = tau[1] - tau[0]
S_XX_fft = np.fft.fftshift(np.abs(np.fft.fft(np.fft.ifftshift(R_XX))) * dtau)
f_fft = np.fft.fftshift(np.fft.fftfreq(len(tau), d=dtau))
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 自己相関関数のプロット
axes[0].plot(tau, R_XX, 'b-', linewidth=2)
axes[0].set_title("Autocorrelation Function $R_{XX}(\\tau)$")
axes[0].set_xlabel("Lag $\\tau$")
axes[0].set_ylabel("$R_{XX}(\\tau)$")
axes[0].axhline(y=0, color='k', linestyle='--', linewidth=0.5)
axes[0].axvline(x=0, color='k', linestyle='--', linewidth=0.5)
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(-5, 5)
# PSDのプロット
axes[1].plot(f, S_XX, 'r-', linewidth=2, label='Analytical')
axes[1].plot(f_fft, S_XX_fft, 'b--', linewidth=1.5, alpha=0.7, label='FFT (numerical)')
axes[1].set_title("Power Spectral Density $S_{XX}(f)$")
axes[1].set_xlabel("Frequency $f$")
axes[1].set_ylabel("$S_{XX}(f)$")
axes[1].legend()
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(-5, 5)
plt.tight_layout()
plt.savefig("autocorrelation_psd.png", dpi=150, bbox_inches="tight")
plt.show()
# 逆フーリエ変換でR_XX(tau)が復元できることを確認
R_XX_recovered = np.fft.fftshift(
np.real(np.fft.ifft(np.fft.ifftshift(S_XX_fft / dtau)))
) * (f_fft[1] - f_fft[0]) * len(f_fft)
print("ウィーナー=ヒンチンの定理の検証:")
print(f" R_XX(0) = {R_XX[len(R_XX)//2]:.4f}")
print(f" ∫S_XX(f)df ≈ {np.trapz(S_XX, f):.4f}")
print(f" 理論値: {sigma2:.4f}")
このコードでは、指数減衰型の自己相関関数 $R_{XX}(\tau) = \sigma^2 e^{-\alpha |\tau|}$ を持つ弱定常過程について、ウィーナー=ヒンチンの定理を数値的に検証しています。対応するPSDは
$$ S_{XX}(f) = \frac{2\sigma^2 \alpha}{\alpha^2 + (2\pi f)^2} $$
という ローレンツ型 の形状になります。この式の導出を行います。
$$ \begin{align} S_{XX}(f) &= \int_{-\infty}^{\infty} \sigma^2 e^{-\alpha|\tau|} e^{-j2\pi f\tau} \, d\tau \\ &= \sigma^2 \left[\int_{-\infty}^{0} e^{\alpha\tau} e^{-j2\pi f\tau} \, d\tau + \int_{0}^{\infty} e^{-\alpha\tau} e^{-j2\pi f\tau} \, d\tau\right] \\ &= \sigma^2 \left[\frac{1}{\alpha – j2\pi f} + \frac{1}{\alpha + j2\pi f}\right] \\ &= \sigma^2 \cdot \frac{(\alpha + j2\pi f) + (\alpha – j2\pi f)}{(\alpha – j2\pi f)(\alpha + j2\pi f)} \\ &= \sigma^2 \cdot \frac{2\alpha}{\alpha^2 + (2\pi f)^2} \end{align} $$
確かにローレンツ型のPSDが得られました。$\alpha$ が小さいほどPSDのピークが鋭くなり、自己相関がゆっくり減衰する（長い記憶を持つ）ことを意味します。
エルゴード性の数値検証
エルゴード性を数値的に確認するシミュレーションを行います。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# AR(1)過程: X[n] = a*X[n-1] + W[n], |a| < 1
# これは弱定常かつ平均エルゴード的
a = 0.9
sigma_w = 1.0
n_total = 100000
# --- サンプルパスを1本生成 ---
X = np.zeros(n_total)
for n in range(1, n_total):
X[n] = a * X[n - 1] + np.random.normal(0, sigma_w)
# 理論値
mu_theory = 0.0
var_theory = sigma_w**2 / (1 - a**2)
R_theory = lambda tau: var_theory * a**np.abs(tau)
# --- 時間平均 vs アンサンブル平均 ---
# 時間平均を様々なTで計算
T_values = np.logspace(1, 4.9, 50).astype(int)
time_means = [np.mean(X[:T]) for T in T_values]
# アンサンブル平均（複数の独立なサンプルパスから）
n_ensemble = 500
ensemble_length = 1000
ensemble_means = []
for _ in range(n_ensemble):
X_ens = np.zeros(ensemble_length)
for n in range(1, ensemble_length):
X_ens[n] = a * X_ens[n - 1] + np.random.normal(0, sigma_w)
ensemble_means.append(np.mean(X_ens[100:])) # バーンイン期間を除く
ensemble_mean = np.mean(ensemble_means)
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (1) 時間平均の収束
axes[0].semilogx(T_values, time_means, 'b-', linewidth=1)
axes[0].axhline(y=mu_theory, color='r', linestyle='--', linewidth=2, label=f'Theory: μ={mu_theory}')
axes[0].axhline(y=ensemble_mean, color='g', linestyle=':', linewidth=2, label=f'Ensemble: {ensemble_mean:.4f}')
axes[0].set_title("Ergodicity: Time Average → Ensemble Average")
axes[0].set_xlabel("Sample size T")
axes[0].set_ylabel("Time average of X")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# (2) 自己相関関数の比較
max_lag = 50
lags = np.arange(0, max_lag + 1)
# 時間平均による推定
R_time = np.array([np.mean(X[:5000] * np.roll(X[:5000], -k))
for k in range(max_lag + 1)])
R_time_normalized = R_time / R_time[0] * var_theory
# 理論値
R_theo_vals = R_theory(lags)
axes[1].plot(lags, R_theo_vals, 'r-', linewidth=2, label='Theory')
axes[1].plot(lags, R_time_normalized, 'b--', linewidth=1.5, label='Time average estimate')
axes[1].set_title("Autocorrelation: Theory vs Time Average")
axes[1].set_xlabel("Lag k")
axes[1].set_ylabel("$R_{XX}(k)$")
axes[1].legend()
axes[1].grid(True, alpha=0.3)
# (3) アンサンブル平均のヒストグラム
axes[2].hist(ensemble_means, bins=30, density=True, alpha=0.7, color='steelblue',
edgecolor='black', linewidth=0.5)
axes[2].axvline(x=mu_theory, color='r', linestyle='--', linewidth=2, label=f'Theory: μ={mu_theory}')
axes[2].set_title(f"Ensemble Mean Distribution (N={n_ensemble} paths)")
axes[2].set_xlabel("Sample mean")
axes[2].set_ylabel("Density")
axes[2].legend()
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("ergodicity_verification.png", dpi=150, bbox_inches="tight")
plt.show()
print(f"AR(1)過程 (a={a}):")
print(f" 理論平均: {mu_theory}")
print(f" 時間平均 (T={n_total}): {np.mean(X):.6f}")
print(f" アンサンブル平均 ({n_ensemble}パス): {ensemble_mean:.6f}")
print(f" 理論分散: {var_theory:.4f}")
print(f" 標本分散 (時間): {np.var(X):.4f}")
まとめ
本記事では、確率過程の定義と分類について包括的に解説しました。
- 確率過程は確率変数の族 $\{X(t), t \in T\}$ として定義され、パラメータ集合 $T$ と状態空間 $S$ による4分類がある
- 定常過程には強定常と弱定常があり、ガウス過程では両者が一致する
- 弱定常過程は平均関数と自己相関関数で2次の統計的性質が完全に特徴づけられる
- エルゴード性は「時間平均 = アンサンブル平均」を保証する重要な性質である
- ウィーナー=ヒンチンの定理により、自己相関関数とパワースペクトル密度がフーリエ変換対をなす
- 主要な確率過程（ベルヌーイ、ランダムウォーク、マルコフ、ポアソン、ガウス、ブラウン運動）は互いに密接に関連している
次のステップとして、以下の記事も参考にしてください。