ポアソン過程

コールセンターにかかってくる電話、Webサイトへのアクセス、放射性物質から飛び出す粒子。いつ起こるかは予測できないが、長い目で見ると一定の割合で起こる出来事は身の回りにあふれている。こうした「ランダムに起こる出来事の回数」を、時間とともに数えていく確率過程がポアソン過程（Poisson process）である。

その名前は、フランスの数学者シメオン・ドニ・ポアソン（Siméon Denis Poisson）が1837年の著作『刑事事件および民事事件における判決の確率に関する研究』で導いたポアソン分布に由来する。ただし、ポアソン自身はこの確率過程を研究したわけではない。19世紀の終わりには、ボルトキェヴィッチ（Ladislaus Bortkiewicz）が1898年の著書『少数の法則』で、プロイセン陸軍の14の騎兵隊について20年間にわたって集計された「馬に蹴られて死亡した兵士の数」がポアソン分布に従うことを示し、ポアソン分布への関心を呼び戻した。

時間の流れの中で起こる出来事の回数をモデル化したのは、デンマークのコペンハーゲン電話会社に勤めていたアーラン（Agner Krarup Erlang）である。アーランは1909年、一定時間内にかかってくる電話の本数がポアソン分布に従うことを示した。この論文は、現在の待ち行列理論につながる最初の論文とされ、電話回線の通信量を表す単位「アーラン」は彼の名にちなんでいる。翌1910年には、ラザフォードとガイガーが放射性物質から放出されるα粒子を数える実験の結果を発表している。なお、「ポアソン過程」という呼び名を最初に文献で用いたのは、1940年のウィリアム・フェラーの論文とされる。

ブラウン運動が連続的に揺れ動く量のモデルだったのに対し、ポアソン過程は「ときどき起こる出来事」を1つずつ数えていく、階段状に増える確率過程である。本記事では、ポアソン過程の定義と基本的な性質、二項分布の極限としての意味、到着間隔が指数分布に従うという見方、そしてデータから強度を推定する方法までを解説する。

ポアソン過程の定義

時刻0から時刻 t までに起こった出来事（イベント）の回数を N_t とする。N_t は 0, 1, 2, \dots の値をとり、イベントが起こるたびに1ずつ増えていく。N_0 = 0 から出発する確率過程 N = (N_t)_{t \geq 0} が次の2つの性質を満たすとき、N を強度 \lambda のポアソン過程という。

独立定常増分は、ブラウン運動の記事で説明したとおり「重ならない区間の増分は互いに独立で、増分の分布は区間の長さだけで決まる」という性質である。イベントの回数でいえば、ある時間帯にイベントが何回起こったかは別の時間帯の回数に影響せず、同じ長さの時間帯なら、いつであっても回数の分布は同じということである。定常増分性と (2) をあわせると、長さ h の任意の区間のイベント数について次が成り立つ。

ポアソン分布 Po(m) の平均と分散はどちらも m であるから

となる。特に \lambda = E[N_1] は単位時間あたりの平均イベント数であり、強度（intensity）や到着率と呼ばれる。たとえば1時間に平均4件の問い合わせが来る窓口なら、時間の単位を「時」として \lambda = 4 である。

図1は \lambda = 1 のポアソン過程のパスの例である。パスは階段状で、イベントが起こった瞬間に1だけ跳ね上がり、それ以外の時間は一定の値にとどまる。跳ね上がる時刻 T_1, T_2, \dots がイベントの発生時刻であり、その間隔 W_1, W_2, \dots はばらばらである。それでも長い時間で見ると、パスは平均 \lambda t の直線のまわりを進んでいく。

白抜きの点は跳ね上がる直前の値、塗りつぶしの点は跳ね上がった後の値を表しており、ちょうど時刻 T_k では N_{T_k} = k となる。

異なる時刻の値の共分散も、ブラウン運動と同じ方法で求められる。s < t のとき N_t = N_s + (N_t - N_s) と分けると、独立増分性より N_s と N_t - N_s は独立なので

となる。

なぜポアソン分布になるのか

定義の (2) でいきなりポアソン分布が出てきたが、これには自然な理由がある。区間 [0, t] を n 個の短い小区間に分けて考えてみよう（図2）。

1つの小区間の長さは \dfrac{t}{n} である。小区間が十分に短ければ、その中でイベントが2回以上起こることはほとんどなく、「1回起こる」か「起こらない」かのどちらかと考えてよい。1回起こる確率は小区間の長さに比例して \lambda \cdot \dfrac{t}{n} とし、各小区間でイベントが起こるかどうかは互いに独立だとする。すると、[0, t] のイベント数 N_t は「成功確率 \dfrac{\lambda t}{n} のベルヌーイ試行を n 回くり返したときの成功回数」となり、二項分布に従う。

小区間を細かくしていく（n \to \infty）と、平均 n \times \dfrac{\lambda t}{n} = \lambda t を一定に保ったまま成功確率が0に近づくので、二項分布はポアソン分布に近づく（ポアソンの少数法則）。

\lambda t = 2 の場合に、分割数 n を増やしたときの確率を比べると次のようになる。n = 100 の時点で、すでにポアソン分布とほとんど変わらない。

つまりポアソン過程は、「短い時間にイベントが起こる確率は時間の長さに比例し、別々の時間帯の出来事は互いに独立」という、ごく自然な仮定から導かれるモデルである。まれに起こる出来事の回数の多くがポアソン分布で近似できるのはこのためである。

到着時刻と到着間隔で見るポアソン過程

ランダムに起こる出来事を記録するには、2通りのやり方がある。たとえばお店の来店者なら、入口のカウンターで「これまでに何人来たか」を数えていく方法と、来店した時刻を「1人目は10時2分、2人目は10時9分、…」と1人ずつ書き留めていく方法である。前者は回数の記録、後者は時刻の記録で、どちらも同じ出来事を別の形で表したものにすぎない。

ここまでのポアソン過程は、回数 N_t の側から見てきた。この節では時刻の側から見直し、2つの見方がどうつながっているかを確かめる。

回数の記録と時刻の記録（計数過程）

k 回目のイベントが起こった時刻を T_k と書く（T_0 = 0 とする）。たとえば受付開始から2分後、9分後、12分後、20分後に問い合わせが来たとすると、T_1 = 2, \; T_2 = 9, \; T_3 = 12, \; T_4 = 20 である。このとき各時刻までの件数 N_t は次のようになる。

たとえば10分までの件数は、到着時刻のリスト 2, 9, 12, 20 のうち10以下のもの（2と9）を数えて N_{10} = 2 となる。つまり時刻の記録さえあれば、どの時刻 t の件数も「t 以下の到着時刻の個数」を数えるだけで求まる。式で書くと

である。ここで I(A) は A が成り立てば1、成り立たなければ0をとる関数（定義関数）で、T_k \leq t を満たす到着時刻だけが1として数えられる。このように、到着時刻の列から件数を数えて作る確率過程を計数過程（counting process）という。「計数」は数を数えるという意味である。到着時刻を時間軸の上に並んだ点の集まりとみて、点過程（point process）と呼ぶこともある。

また、前のイベントから次のイベントまでの待ち時間

を到着間隔という。上の例では W_1 = 2, \; W_2 = 7, \; W_3 = 3, \; W_4 = 8 である。逆に、到着時刻は待ち時間を順に足したもので、T_k = W_1 + W_2 + \cdots + W_k となる。図1の下の帯に描いた T_1, T_2, T_3 と W_1, W_2, W_3 は、この関係を表している。

到着間隔は指数分布に従う

では、ポアソン過程の到着間隔はどんな分布に従うのだろうか。まず最初の到着までの待ち時間 W_1 = T_1 を考える。「最初の到着が時刻 t より後」ということは「時刻 t までに1回も起きていない」ということなので、ポアソン分布の k = 0 の確率から

となる。累積分布関数で書くと P(W_1 \leq t) = 1 - e^{-\lambda t} であり、これは平均 \dfrac{1}{\lambda} の指数分布 Exp(\lambda) の累積分布関数そのものである。つまり、最初の到着までの待ち時間は指数分布に従う。

2回目以降の待ち時間も同じである。独立定常増分性により、イベントが起きた直後から先の回数の増え方は、それまでの経緯とは独立で、時刻0から見たときと同じように振る舞う。いわばイベントが起きるたびに過程がリセットされるので、W_2, W_3, \dots も W_1 と同じ Exp(\lambda) に従い、互いに独立になることが知られている。

すると、k 回目の到着時刻 T_k = W_1 + \cdots + W_k は、独立な指数分布の和になる。指数分布 Exp(\lambda) は形状母数1のガンマ分布 Ga\left(1, \; \dfrac{1}{\lambda}\right) であり、ガンマ分布には「独立な確率変数の和では形状母数が足し算になる」という再生性があるので

となる。平均は E[T_k] = \dfrac{k}{\lambda} で、平均 \dfrac{1}{\lambda} の待ち時間を k 回分足したものになっている。図3は \lambda = 1 のときの T_1, T_2, T_3 の確率密度関数である。

T_1 は指数分布で、待ち時間0の付近が最も起こりやすい。k が大きくなるほど分布は右にずれて広がっていく。

逆に、到着間隔から回数の分布を求める

今度は逆向きに、到着間隔 W_1, W_2, \dots が互いに独立に Exp(\lambda) に従うなら、回数 N_t はポアソン分布 Po(\lambda t) に従うことを確かめよう。これが示せれば、ポアソン過程は「回数がポアソン分布に従う過程」と見ても、「到着間隔が指数分布に従う過程」と見てもよいことになる。

手順1：「ちょうど k 回」を到着時刻で言い換える。 時刻 t までにちょうど k 回起こるのは、k 回目が時刻 t までに起こり（T_k \leq t）、そこから次の到着までの待ち時間 W_{k+1} が、残り時間 t - T_k より長いときである。

手順2：k = 1 で計算してみる。 最初の到着が時刻 s（0 \leq s \leq t）に起こり、残りの t - s の間に次の到着が起こらない場合を考える（図4）。

T_1 = W_1 の密度は \lambda e^{-\lambda s}、W_2 が t - s より長い確率は e^{-\lambda(t-s)} である。両者は独立なので掛け合わせ、最初の到着時刻 s としてありうるすべての値について積分する。

ポイントは、e^{-\lambda s} \cdot e^{-\lambda (t - s)} = e^{-\lambda t} と、指数の部分が最初の到着時刻 s によらない形にまとまることである。そのため積分は「s のとりうる範囲の長さ t」をかけるだけになり、ポアソン分布の k = 1 の確率が現れた。

手順3：一般の k で計算する。 k 回目の到着時刻 T_k の密度 f_k(s) を使えば、k = 1 とまったく同じ計算ができる。W_{k+1} は T_k = W_1 + \cdots + W_k と独立なので

となり、N_t \sim Po(\lambda t) が確かめられた（最後に (k-1)! \times k = k! を使った）。ここでも e^{-\lambda s} \cdot e^{-\lambda(t-s)} = e^{-\lambda t} とまとまるおかげで、残る積分は s^{k-1} の積分だけになっている。なお k = 0 の場合は、前の項で求めた P(N_t = 0) = P(W_1 > t) = e^{-\lambda t} がそのまま使える。

計算例：問い合わせ窓口

30分間に1件も来ない確率

30分間（t = 0.5）の件数は N_{0.5} \sim Po(4 \times 0.5) = Po(2) なので

となる。到着間隔で考えても、最初の問い合わせまでの時間 W_1 \sim Exp(4) が0.5時間を超える確率なので P(W_1 > 0.5) = e^{-4 \times 0.5} = e^{-2} となり、同じ答えが得られる。

1時間にちょうど4件来る確率

平均どおりの4件ちょうどになる確率は、意外にも2割に満たない。

3件目が30分以内に来る確率

「3件目の到着時刻 T_3 が0.5時間以内」は「30分間に3件以上来る」と同じことなので

となる。ガンマ分布 T_3 \sim Ga\left(3, \; \dfrac{1}{4}\right) の累積分布関数を直接計算しても同じ値になる。このように「k 回目の到着が t 以前」と「時刻 t までに k 回以上起こる」は同じ事象であり、どちらで計算しやすいかで使い分けるとよい。また、3件目までの平均待ち時間は E[T_3] = \dfrac{3}{\lambda} = 0.75 時間、すなわち45分である。

たくさん来た後でも次の来やすさは変わらない

最初の30分で3件の問い合わせがあったとする。次の30分に1件も来ない確率は、独立増分性により最初の30分の結果に影響されないので、やはり e^{-2} \approx 0.135 のままである。「たくさん来たから、しばらくは来ないだろう」という直感は、ポアソン過程では成り立たない。

パラメータ推定

強度 \lambda が未知のとき、データから推定することを考えよう。データの取り方として、次の2つの形式がある。

(I) イベントの発生時刻 T_1, T_2, \dots, T_n を記録する。
(II) 一定の時間間隔 \Delta ごとに、その間のイベント数を記録する。

(I) ではイベントの正確な時刻がわかるが、(II) では回数しかわからない。また、(I) の n はイベントの回数、(II) の n は観測区間の個数であることに注意しよう。どちらの場合も、最尤推定法の手順どおり、対数尤度関数を書いて微分し、0とおいて解く。

(I) 発生時刻を記録した場合

到着間隔 W_k = T_k - T_{k-1} に変換すると、W_1, \dots, W_n は互いに独立に指数分布 Exp(\lambda)（密度 f(w) = \lambda e^{-\lambda w}）に従う。尤度関数は各密度の積なので

である。到着間隔の和は最後の到着時刻 \displaystyle\sum_{k=1}^{n} W_k = T_n になるので、対数をとると

となる。これを \lambda で微分して0とおくと

を得る。\dfrac{d\ell}{d\lambda} は \lambda < \hat{\lambda} で正、\lambda > \hat{\lambda} で負なので、\hat{\lambda} で最大となる。また、モーメント法で E[W_k] = \dfrac{1}{\lambda} を標本平均で置き換えても、\dfrac{1}{\hat{\lambda}} = \dfrac{1}{n}\displaystyle\sum_{k=1}^{n} W_k = \dfrac{T_n}{n} より同じ推定量が得られる。

(II) 一定間隔で件数を記録した場合

時間間隔 \Delta ごとの件数を

とおくと、独立定常増分性より M_1, \dots, M_n は互いに独立に Po(\lambda\Delta) に従う。確率関数 P(M_k = m) = e^{-\lambda\Delta}\dfrac{(\lambda\Delta)^m}{m!} の積の対数をとると

となる。\log(\lambda\Delta) = \log \lambda + \log \Delta に注意して \lambda で微分し、0とおくと

を得る。最後の等号では、件数の和が途中の項で打ち消し合って \displaystyle\sum_{k=1}^{n} M_k = N_{n\Delta} - N_0 = N_{n\Delta} となることを使った。

計算例

n = 8, \; T_8 = 40 より

すなわち1時間あたり12件と推定される。このデータの対数尤度関数 \ell(\lambda) = 8 \log \lambda - 40\lambda を描くと図5のようになり、\lambda = 0.2 で山の頂上（接線の傾きが0）になっていることが確かめられる。

合計は 7 + 4 + 6 + 9 + 5 + 3 + 8 + 6 = 48 件、観測時間は 5 \times 8 = 40 分なので

と推定される。

練習問題

まとめ

Python実装

到着間隔が指数分布に従うという性質を使ってポアソン過程を生成し、確率の計算、強度の推定、性質の確認を行う。

import numpy as np
from scipy.stats import poisson, gamma

def simulate_poisson_process(lam, T, rng):
    """強度 lam のポアソン過程の到着時刻を区間 [0, T] で生成する"""
    times = []
    t = rng.exponential(1 / lam)          # 到着間隔は平均 1/lam の指数分布
    while t <= T:
        times.append(t)
        t += rng.exponential(1 / lam)
    return np.array(times)

def estimate_from_times(times):
    """(I) 到着時刻 T_1, ..., T_n から λ の最尤推定値 n / T_n を求める"""
    return len(times) / times[-1]

def estimate_from_counts(counts, dt):
    """(II) 間隔 dt ごとの件数 M_1, ..., M_n から λ の最尤推定値を求める"""
    return np.sum(counts) / (len(counts) * dt)

# ========== 計算例：問い合わせ窓口（λ = 4 件/時） ==========
print("【問い合わせ窓口（λ = 4 件/時）】")
lam = 4.0
print(f"30分間に1件も来ない確率: {poisson.pmf(0, lam * 0.5):.4f}")
print(f"1時間にちょうど4件来る確率: {poisson.pmf(4, lam * 1):.4f}")
print(f"3件目が30分以内に来る確率: {poisson.sf(2, lam * 0.5):.4f}")    # P(N ≥ 3)
print(f"  （ガンマ分布 T_3 で計算しても同じ）: {gamma.cdf(0.5, a=3, scale=1 / lam):.4f}")

# ========== パラメータ推定の例 ==========
print("\n【パラメータ推定】")
times = np.array([2, 9, 12, 20, 21, 27, 35, 40])        # 到着時刻（分）
print(f"(I) 到着時刻から: λ の推定値 = {estimate_from_times(times):.4f} 件/分")
counts = np.array([7, 4, 6, 9, 5, 3, 8, 6])             # 5分ごとの件数
print(f"(II) 5分ごとの件数から: λ の推定値 = {estimate_from_counts(counts, 5):.4f} 件/分")

# ========== シミュレーションで性質を確認 ==========
print("\n【シミュレーション（λ = 2, T = 1000）】")
rng = np.random.default_rng(1)
times = simulate_poisson_process(2.0, 1000, rng)
gaps = np.diff(np.concatenate([[0], times]))
print(f"イベント数: {len(times)}, 到着間隔の平均: {gaps.mean():.4f}（理論値 1/λ = 0.5）")
counts = np.histogram(times, bins=np.arange(0, 1001, 5))[0]   # 長さ5の区間ごとの件数
print(f"長さ5の区間の件数: 平均 {counts.mean():.3f}, 分散 {counts.var():.3f}（理論値はどちらも λt = 10）")
print(f"(I) の推定値: {estimate_from_times(times):.4f}, (II) の推定値: {estimate_from_counts(counts, 5):.4f}")

# ========== (I) の推定量の偏り ==========
print("\n【(I) の推定量 n/T_n の期待値（λ = 0.2, n = 8, 10万回）】")
T8 = rng.gamma(shape=8, scale=1 / 0.2, size=100000)     # 8件目の到着時刻 T_8
print(f"n/T_n の平均: {np.mean(8 / T8):.4f}（理論値 nλ/(n-1) = {8 * 0.2 / 7:.4f}）")
print(f"(n-1)/T_n の平均: {np.mean(7 / T8):.4f}")

シミュレーションでは、到着間隔の平均が理論値 \dfrac{1}{\lambda} = 0.5 に近く、長さ5の区間の件数の平均と分散はどちらも理論値 \lambda t = 10 に近い値になっている。平均と分散がほぼ等しいことは、件数がポアソン分布に従っていることを示している。また、(I) の推定量 \dfrac{n}{T_n} を n = 8 で10万回計算した平均は理論値 \dfrac{n}{n-1}\lambda \approx 0.2286 とほぼ一致し、\dfrac{n-1}{T_n} にすると真の値0.2にほぼ一致する。