ブラウン運動

1827年の夏、スコットランドの植物学者ロバート・ブラウン（Robert Brown）は、水に浸した植物（Clarkia pulchella）の花粉を顕微鏡で観察していて、花粉に含まれる微粒子が休むことなく不規則に動き回っていることに気づいた。まるで生きているかのような動きだったが、石炭の粉やガラス、金属の粉末でも同じ動きが見られたことから、生命とは関係のない現象であることがわかった。この観察は翌1828年に発表され、のちにこの不規則な運動は彼の名をとってブラウン運動（Brownian motion）と呼ばれるようになった。

この不規則な動きを数式で表す試みは、意外にも物理学ではなく金融の世界から始まった。フランスのルイ・バシュリエ（Louis Bachelier）は1900年、アンリ・ポアンカレの指導のもとで書いた博士論文『投機の理論』（Théorie de la spéculation）で、株式市場の価格変動を現在のブラウン運動にあたる確率モデルで表し、オプションの価格評価に用いた。その5年後の1905年、アルベルト・アインシュタインはブラウン運動を、水の分子が微粒子に絶えず衝突する結果として説明し、微粒子の移動距離が経過時間ではなく経過時間の平方根に比例して広がることを理論的に導いた。この予測は1908〜1909年の実験で確かめられ、原子や分子が実在することを裏づける証拠となった（実験を行ったジャン・ペランは、物質の不連続な構造に関する研究で1926年にノーベル物理学賞を受賞している）。

ブラウン運動に初めて完全で厳密な数学的解析を与えたのは、アメリカの数学者ノーバート・ウィーナー（Norbert Wiener）で、1923年のことである。このため数学ではブラウン運動をウィーナー過程（Wiener process）とも呼ぶ。現在では、物理学の拡散現象から株価や為替レートのモデルまで、時間とともにランダムに揺れ動く量を表す最も基本的な確率モデルとして広く使われている。

本記事では、マルコフ連鎖に続く確率過程の基本として、確率過程と独立定常増分の考え方から始め、ブラウン運動の定義と性質、中心極限定理を通じたランダムウォークとの関係、そして一定の間隔で観測したデータからパラメータを推定する方法までを解説する。

確率過程とパス

1つの確率変数は「1回の観測でどんな値が出るか」のばらつきを表す。これに対して、微粒子の位置や株価のように時間とともに変化していくランダムな量を表すには、時刻ごとに確率変数を用意し、それらを時刻の順に並べたものを考えればよい。

各時刻 t \geq 0 に対して確率変数 X_t が与えられているとき、その集まり

を確率過程（stochastic process）という。添え字 t は時間を表すことが多い。t が実数のように連続的な値をとる場合を連続時間確率過程、t = 0, 1, 2, \dots のように飛び飛びの値をとる場合を離散時間確率過程という。マルコフ連鎖は離散時間の確率過程の代表例であり、本記事で扱うブラウン運動は連続時間の確率過程である。

確率過程を1回実現させると、各時刻の値 x_t が決まり、t \mapsto x_t という時間の関数のグラフが1本描かれる。これを確率過程のパス（path、見本路）という。同じ確率過程でも、実現させるたびに異なるパスが描かれる点が重要である。図1は、後で定義する標準ブラウン運動を3回実現させたパスである。

どのパスも出発点は0だが、その後の動きはまったく異なる。確率過程を学ぶとは、このように無数に描かれうるパス全体がもつ統計的な性質、たとえば各時刻の値の分布や、異なる時刻の値どうしの関係を調べることにほかならない。

独立定常増分

確率過程の性質を調べるときに中心となるのが、ある時刻から別の時刻までの変化量

である。これを増分（increment）という。応用上重要な確率過程の多くは、増分について次の2つの性質をもつ。

(1) は、重なり合わない時間区間での変化量どうしが互いに影響しないことを意味する。前の区間で大きく上がったからといって、次の区間で下がりやすくなったり上がりやすくなったりすることはない。(2) は、変化量の分布が区間の長さ h だけで決まり、どの時刻から測り始めたかによらないことを意味する。この2つを満たす確率過程を独立定常増分過程（process of independent and stationary increments）という。ブラウン運動のほか、ランダムに起こるイベントの回数を数えるポアソン過程もこの性質をもつ。

ブラウン運動の定義

B_0 = 0 から出発する確率過程 B = (B_t)_{t \geq 0} が次の3つの性質を満たすとき、B をブラウン運動という。

(2) は、時刻 t における位置が正規分布に従い、その平均も分散も経過時間 t に比例して大きくなることを表している。パラメータ \mu は平均的に進む向きと速さを表し、ドリフト（drift）と呼ばれる。\sigma > 0 は揺れの大きさを表し、金融の分野ではボラティリティ（volatility）と呼ばれる。(3) は、パスが途切れたり飛び移ったりしないことを表す。微粒子が瞬間移動することはないので、物理的にも自然な条件である。

定義 (1) の定常増分性から、s < t のとき増分 B_t - B_s は B_{t-s} - B_0 = B_{t-s} と同じ分布に従う。これと (2) をあわせると、次の結果が得られる。

すなわち、長さ h の時間に生じる変化は平均 \mu h、分散 \sigma^2 h の正規分布に従い、しかも重ならない区間どうしの変化は独立である。ブラウン運動に関する計算のほとんどは、この性質だけで行うことができる。

標準ブラウン運動（ウィーナー過程）

特に \mu = 0, \; \sigma^2 = 1 の場合、すなわち W_t \sim N(0, t) となるものを標準ブラウン運動またはウィーナー過程（ウィナー過程とも表記される）といい、W = (W_t)_{t \geq 0} と書くことが多い。図1は標準ブラウン運動のパスである。

一般のブラウン運動は、標準ブラウン運動 W を使って次のように表すことができる。

右辺の \mu t は「一定の速さで進む直線的な動き」、\sigma W_t は「その周りのランダムな揺れ」である。実際、W の増分 W_{t+h} - W_t \sim N(0, h) を使うと

となり、増分の分布は区間の長さ h だけで決まる（定常増分性）。独立増分性とパスの連続性も W からそのまま引き継がれるので、B_t = \mu t + \sigma W_t はブラウン運動の定義をすべて満たす。

ばらつきは時間の平方根に比例して広がる

B_t の平均は \mu t、分散は \sigma^2 t なので、標準偏差は \sigma \sqrt{t} である。正規分布の性質から、各時刻 t において

が成り立つ。図2は \mu = 0.5, \; \sigma = 1 のブラウン運動について、平均 \mu t とこの95%の範囲を描き、パスの例を重ねたものである。

平均は直線的に増えていくのに対し、95%の範囲は放物線を横に倒したような形で広がっていく。たとえば t = 1 では範囲は [-1.46, \; 2.46] だが、t = 4 では [-1.92, \; 5.92]、t = 9 では [-1.38, \; 10.38] となり、範囲の幅は時間が4倍で2倍、9倍で3倍になる。なお、この範囲は各時刻ごとに値が入る確率が95%という意味であり、パス全体がこの範囲に収まる確率が95%という意味ではない。

異なる時刻の値の共分散

ブラウン運動の時刻 s と t（s < t）の値の共分散を求めてみよう。ポイントは、B_t を「時刻 s までの値」と「その後の増分」に分けることである。

独立増分性より B_s = B_s - B_0 と B_t - B_s は独立なので、その共分散は0である。したがって

となる。s と t の大小が逆の場合も同様なので、まとめて次のように書ける。

共分散が正になるのは、B_t が B_s を「出発点」として含んでいるからである。たとえば標準ブラウン運動では \mathrm{Cov}[W_2, W_5] = 2 であり、相関係数は \dfrac{2}{\sqrt{2 \times 5}} \approx 0.632 となる。一般に標準ブラウン運動の相関係数は \sqrt{\dfrac{s}{t}} であり、2つの時刻が近いほど値は強く相関し、離れるほど相関は弱くなる。

ランダムウォークとの関係

ブラウン運動は連続時間の確率過程だが、時間を細かく区切って考えると理解しやすい。区間 [0, 1] を n 等分して t_k = \dfrac{k}{n}（k = 0, 1, \dots, n）とおくと、

と書ける。定義より各増分 \varepsilon_i は互いに独立に N\left(\dfrac{\mu}{n}, \; \dfrac{\sigma^2}{n}\right) に従う。つまり、ブラウン運動を等間隔に観測した値は、独立な確率変数を次々に足していった和、すなわちランダムウォーク（random walk）になっている。分割数 n を大きくするほど、点を結んだ折れ線はブラウン運動のパスに近づいていく。

興味深いのは、1歩ごとの増分が正規分布でなくても同じことが起こる点である。先ほどのコイン投げゲームで、1歩の大きさを \pm\dfrac{1}{\sqrt{n}} に縮めて、時間1の間に n 歩進む場合を考える。1歩の平均は0、分散は \dfrac{1}{n} なので、時刻1での位置（n 歩の合計）は平均0、分散1となり、中心極限定理により n \to \infty で標準正規分布 N(0, 1) に近づく。図3は n = 10, \; 100, \; 1000 の場合のパスである。

n = 10 ではカクカクした折れ線にすぎないが、n = 1000 では図1の標準ブラウン運動とほとんど見分けがつかない。各時刻の値だけでなく、パス全体としても標準ブラウン運動に近づくことが知られており、これをドンスカーの不変原理（関数中心極限定理）という。ブラウン運動が自然界や社会のさまざまな場面に現れるのは、小さなランダムな変化が独立に積み重なる現象であれば、個々の変化の分布の細かい形によらず、その累積が全体としてブラウン運動で近似できるからである。

計算例：試合の得点差

ブラウン運動の独立定常増分性を使うと、「途中経過がわかったときに最終結果を予測する」計算が簡単にできる。実際に、スポーツの試合の得点差の推移をブラウン運動でモデル化した研究もある（H. S. Stern, 1994, Journal of the American Statistical Association）。ここでは数値を単純にした例で考えよう。

試合開始前にホームチームが勝つ確率

試合終了時の得点差は B_1 \sim N(2, \; 10^2) である。ホームチームが勝つのは B_1 > 0 のときなので、標準化すると

となる。ここで Z は標準正規分布に従う確率変数、\Phi はその累積分布関数である。試合前の時点では、ホームチームが勝つ確率は約58%にすぎない。

第3クォーター終了時に4点リードしている場合

第3クォーター終了時（t = 0.75）に、ホームチームが4点リードしている（B_{0.75} = 4）とする。最終的な得点差を「現在のリード」と「残り時間に生じる増分」に分けると

となる。独立増分性により、残り時間の増分はそれまでの試合展開とは独立であり、定常増分性によりその分布は残り時間の長さ0.25だけで決まる。

よって最終的な得点差は N(4.5, \; 5^2) に従い、ホームチームが勝つ確率は

となる。

同じ4点リードでも時刻によって意味が違う

ハーフタイム（t = 0.5）に4点リードしている場合は、残り時間が0.5なので増分は N(1, \; 10^2 \times 0.5) = N(1, \; 50) に従い、最終的な得点差は N(5, \; 50) に従う。\sqrt{50} \approx 7.07 より

である。同じ4点のリードでも、残り時間が短いほど逆転に必要な大きな揺れが起こりにくいため、勝つ確率は高くなる。残り時間に生じる揺れの標準偏差が \sigma\sqrt{1 - t} と、残り時間の平方根に比例して小さくなることがこの差を生んでいる。

パラメータ推定

ブラウン運動のパスは時間について連続だが、実際のデータは「1秒ごと」「1日ごと」のように飛び飛びの時刻でしか観測できない。そこで、ブラウン運動 B_t \sim N(\mu t, \; \sigma^2 t) を時間間隔 \Delta > 0 で観測して

というデータを得たとき、未知のパラメータ \mu, \; \sigma^2 を推定する方法を考えよう。観測期間の長さを T = n\Delta とおく。

増分に変換する

観測値 B_{k\Delta} どうしは互いに相関しているので、そのままでは扱いにくい。そこで隣り合う観測値の差（増分）

に変換する。独立定常増分性より、Z_1, \dots, Z_n は互いに独立に同じ正規分布 N(\mu\Delta, \; \sigma^2\Delta) に従う。これで問題は「正規分布からの無作為標本にもとづく推定」という、おなじみの形に帰着した。

最尤推定法では、観測されたデータが得られる確率（密度）をパラメータの関数とみて尤度関数と呼び、それを最大にするパラメータの値を推定値とする。最大となる点は、関数を各パラメータで偏微分して0とおいた方程式を解いて求める。以下、この計算を3つの手順に分けて進めよう。

手順1：対数尤度関数を書く

1つの増分 Z_k は平均 \mu\Delta、分散 \sigma^2\Delta の正規分布に従うので、正規分布の確率密度関数の平均と分散のところに \mu\Delta と \sigma^2\Delta を入れて

となる。Z_1, \dots, Z_n は独立なので、データ全体の同時確率密度は各密度の積であり、これが尤度関数である。

積のままでは微分しにくいので、対数をとって和に直す。対数は単調増加関数なので、L を最大にするパラメータと \log L を最大にするパラメータは同じである。\log \dfrac{1}{\sqrt{a}} = -\dfrac{1}{2}\log a と \log e^{x} = x を使うと、1つの増分あたり

となる。これを k = 1, \dots, n について足し合わせると、第1項は k によらないので n 倍になり、次の対数尤度関数が得られる。

第1項は \sigma^2 だけを含む。第2項は「各増分 Z_k と、その期待値 \mu\Delta とのずれの二乗和」に負の係数をかけたものである。

手順2：μ で偏微分して0とおく

\mu を含むのは第2項だけである。合成関数の微分 \dfrac{\partial}{\partial \mu}(Z_k - \mu\Delta)^2 = 2(Z_k - \mu\Delta) \times (-\Delta) を使うと

となる。これを0とおくと \displaystyle\sum_{k=1}^{n} Z_k = n\hat{\mu}\Delta より

を得る。ここで \bar{Z} は増分の平均である。偏微分は \mu < \hat{\mu} のとき正、\mu > \hat{\mu} のとき負になるので、\ell は \hat{\mu} まで増えてその後は減る。つまり \hat{\mu} で最大となる。また、この答えには \sigma^2 が含まれていない。

さらに、増分の和は途中の項が打ち消し合って

となるので、\hat{\mu} = \dfrac{B_{n\Delta} - B_0}{T} とも書ける。ドリフトの推定値は「全体の変化量 ÷ 観測期間」であり、途中の観測値にはまったく依存しない。

手順3：σ² で偏微分して0とおく

手順2の答えは \sigma^2 によらないので、\sigma^2 がどんな値でも \mu = \hat{\mu} のときに \ell は最大になる。そこで \mu = \hat{\mu}（つまり \mu\Delta = \bar{Z}）を代入してから、\sigma^2 について最大化すればよい。このとき二乗和は

となる。\sigma^2 を1つの変数 v とみて、\log(2\pi v\Delta) = \log(2\pi\Delta) + \log v と分けておくと

である。第1項は定数なので、\dfrac{d}{dv} \log v = \dfrac{1}{v}、\dfrac{d}{dv} \dfrac{1}{v} = -\dfrac{1}{v^2} を使って v で微分すると

となる。これを0とおくと \dfrac{S}{\Delta} = nv より

を得る。括弧の中の \dfrac{S}{\Delta} - nv は、v が \hat{\sigma}^2 より小さいと正、大きいと負になる。したがって \ell は \hat{\sigma}^2 まで増えてその後は減るので、ここで最大となる。

以上をまとめると、次の推定量が得られる。

両辺に \Delta をかけると、\hat{\mu}\Delta = \bar{Z} は増分の標本平均、\hat{\sigma}^2\Delta = \dfrac{1}{n}\displaystyle\sum_{k=1}^{n}(Z_k - \bar{Z})^2 は増分の標本分散である。つまりこの推定量は、増分の平均と分散を1区間の長さ \Delta で割って「単位時間あたり」に直したものになっている。

モーメント法でも同じ推定量になる

モーメント法で求めても同じ推定量が得られる。Z_k \sim N(\mu\Delta, \; \sigma^2\Delta) より

であるから、左辺を標本モーメントで置き換えた連立方程式

を解けばよい。第2式から \hat{\sigma}^2\Delta = \dfrac{1}{n}\displaystyle\sum_{k=1}^{n} Z_k^2 - \bar{Z}^2 = \dfrac{1}{n}\displaystyle\sum_{k=1}^{n} (Z_k - \bar{Z})^2 となり、最尤推定量と一致する。

計算例：微粒子の位置データ

観測間隔は \Delta = 2、増分の個数は n = 8、観測期間は T = 16 である。まずドリフトは

と求まる。次に、増分の平均は \bar{Z} = \dfrac{8}{8} = 1 で、平均からの偏差 Z_k - \bar{Z} は 1, \; -2, \; 2, \; -1, \; 0, \; 1, \; -2, \; 1 となるので

となる。微粒子は流れによって平均して毎秒 0.5 μm 進み、それに加えて1秒あたり分散 1.0 μm² の割合でランダムに揺れている、と推定される。

この推定値で対数尤度関数が本当に最大になっているかを、グラフで確かめてみよう。このデータでは n = 8, \; \Delta = 2 なので、対数尤度関数は

である。図4の左は \sigma^2 = 1.0 に固定して \mu を動かしたもの、右は \mu = 0.5 に固定して \sigma^2 を動かしたものである。

左のグラフは \mu の2次関数（上に凸の放物線）で、頂上は \mu = 0.5 にある。右のグラフは \sigma^2 が小さいところで急に下がり、\sigma^2 = 1.0 で最大になったあと、ゆるやかに下がっていく。手順2と手順3で「偏微分を0とおいた」のは、この山の頂上、つまり接線の傾きが0になる点を探していたことにあたる。

推定の精度と観測間隔

2つの推定量の精度は、観測のしかたによって大きく異なる。\hat{\mu} = \dfrac{B_T - B_0}{T} の分散は、B_T - B_0 \sim N(\mu T, \; \sigma^2 T) より

となり、観測期間 T だけで決まって、観測間隔 \Delta にはよらない。一方 \hat{\sigma}^2 については、\dfrac{n\hat{\sigma}^2}{\sigma^2} = \dfrac{1}{\sigma^2\Delta}\displaystyle\sum_{k=1}^{n} (Z_k - \bar{Z})^2 が自由度 n - 1 のカイ二乗分布に従うことから

となり、観測回数 n を増やすほど精度が上がる。観測期間を変えずに観測間隔を細かくすると、\hat{\sigma}^2 はいくらでも正確になるが、\hat{\mu} の精度はまったく改善しない。ドリフトを正確に知るには、より長い期間にわたって観測を続けるしかない。

練習問題

まとめ

Python実装

ブラウン運動のパスの生成と、等間隔の観測データからのパラメータ推定をPythonで実装する。パスの生成では、独立定常増分性にしたがって正規乱数で増分を作り、その累積和をとっている。

import numpy as np
from scipy.stats import norm

def simulate_bm(mu, sigma, T, n, rng):
    """ドリフト mu、拡散 sigma のブラウン運動を区間 [0, T] で n 分割して生成する"""
    dt = T / n
    # 独立定常増分：各区間の増分は N(mu*dt, sigma^2*dt) に独立に従う
    increments = rng.normal(mu * dt, sigma * np.sqrt(dt), size=n)
    t = np.linspace(0, T, n + 1)
    B = np.concatenate([[0.0], np.cumsum(increments)])
    return t, B

def estimate_bm(B, dt):
    """間隔 dt で観測した B_0, B_dt, ..., B_ndt から (mu, sigma^2) の最尤推定値を求める"""
    Z = np.diff(B)                                  # 増分 Z_1, ..., Z_n
    mu_hat = Z.mean() / dt                          # = (B_T - B_0) / T
    sigma2_hat = np.mean((Z - Z.mean())**2) / dt
    return mu_hat, sigma2_hat

# ========== 計算例：試合の得点差 ==========
print("【試合の得点差（μ = 2, σ = 10）】")
mu, sigma = 2.0, 10.0
print(f"試合前にホームチームが勝つ確率: {norm.cdf(mu / sigma):.4f}")
for t, lead in [(0.5, 4), (0.75, 4)]:
    m = lead + mu * (1 - t)          # 最終得点差の平均
    s = sigma * np.sqrt(1 - t)       # 最終得点差の標準偏差
    print(f"t = {t} で {lead} 点リード → 勝つ確率: {norm.cdf(m / s):.4f}")

# ========== 計算例：微粒子の位置データ ==========
print("\n【微粒子の位置データ（Δ = 2 秒）】")
B = np.array([0, 2, 1, 4, 4, 5, 7, 6, 8], dtype=float)
mu_hat, sigma2_hat = estimate_bm(B, dt=2.0)
print(f"μ の推定値: {mu_hat:.4f}, σ² の推定値: {sigma2_hat:.4f}")

# ========== シミュレーションで推定を確認 ==========
print("\n【シミュレーション（真の値 μ = 0.5, σ² = 1）】")
rng = np.random.default_rng(1)
t, B = simulate_bm(mu=0.5, sigma=1.0, T=100, n=10000, rng=rng)
mu_hat, sigma2_hat = estimate_bm(B, dt=t[1] - t[0])
print(f"T = 100, Δ = 0.01: μ の推定値 = {mu_hat:.4f}, σ² の推定値 = {sigma2_hat:.4f}")

# ========== 観測間隔と推定精度 ==========
print("\n【観測期間 T = 10 のまま観測間隔 Δ を変える（各2000回）】")
print(f"理論値: μ の推定値の標準偏差 σ/√T = {1 / np.sqrt(10):.4f}")
for dt in [1.0, 0.1, 0.01]:
    n = int(round(10 / dt))
    est = np.array([estimate_bm(simulate_bm(0.5, 1.0, 10, n, rng)[1], dt)
                    for _ in range(2000)])
    print(f"Δ = {dt:4.2f}, n = {n:4d}: μ の推定値の標準偏差 = {est[:, 0].std():.4f}, "
          f"σ² の推定値の標準偏差 = {est[:, 1].std():.4f}")

シミュレーションでは、\sigma^2 の推定値は真の値1にきわめて近いのに対し、\mu の推定値は0.39と、真の値0.5から0.11ほどずれている。観測期間が T = 100 のとき \hat{\mu} の標準偏差は \dfrac{\sigma}{\sqrt{T}} = 0.1 なので、この程度のずれは珍しくない。最後の実験では、観測期間を固定して観測間隔を細かくしても \hat{\mu} の標準偏差は理論値0.3162の付近から変わらない。一方 \hat{\sigma}^2 の標準偏差は、理論値 \dfrac{\sqrt{2(n-1)}}{n} \approx 0.424, \; 0.141, \; 0.045 とほぼ同じ割合で小さくなっており、「推定の精度と観測間隔」で述べた性質が確認できる。