複合ポアソン過程

保険会社には、毎日のように保険金の請求が届く。何件の請求が来るかは日によって違い、1件あたりの金額も請求ごとに違う。保険会社にとって大切なのは、一定の期間に支払う保険金の総額がどれくらいになり、どれくらいばらつくかである。このように、ランダムに起こる出来事のそれぞれにランダムな大きさがついていて、その合計を時間とともに追いかける確率過程を複合ポアソン過程（compound Poisson process）という。

複合ポアソン過程は保険数学の中から生まれた。スウェーデンのアクチュアリー（保険数理の専門家）フィリップ・ルンドベリ（Filip Lundberg）は、1903年の博士論文で、保険会社が支払う保険金の総額を、現在では複合ポアソン過程と呼ばれるモデルで表した。これはそれまで定義されていなかった新しい概念であり、ルンドベリはこのモデルについての中心極限定理も独自の方法で導いている。ルンドベリの考えは、クラメール・ラオの不等式に名を残すスウェーデンの数学者ハラルド・クラメール（Harald Cramér）とその弟子たちの仕事を通じて広く知られるようになった。クラメールは1930年と1955年の著書で、この理論の発展を明快にまとめている。保険料の収入と保険金の支払いを組み合わせた「クラメール・ルンドベリ・モデル」は、保険会社が破産する確率を調べる破産理論（ruin theory）の基礎となるモデルである。

ポアソン過程は、出来事が起こるたびに1ずつ増える確率過程だった。複合ポアソン過程は、その「1ずつ」を「出来事ごとにランダムな大きさ」に置き換えたものである。保険金の総額のほか、店の売上の合計（来店する客の数も、1人あたりの購入額もランダム）など、「いつ起こるかわからない出来事の大きさの合計」を表すモデルとして使われる。本記事では、複合ポアソン過程の定義、条件付き期待値を使った平均と分散の求め方、モーメント母関数、間引きと呼ばれる重要な特別な場合、正規近似を使った計算例、そしてデータからのパラメータ推定までを解説する。

複合ポアソン過程の定義

具体例から始めよう。あるパン屋には、客が1時間あたり平均8人の割合でランダムに来店する。1人あたりの購入額は客によってばらばらで、平均は1.5千円（1500円）である。時刻 t までに来店した客の数を N_t とし、N = (N_t)_{t \geq 0} は強度 \lambda = 8（人/時）のポアソン過程に従うとする。k 番目に来店した客の購入額を U_k とすると、時刻 t までの売上の合計は

と表せる。たとえば時刻 t までに3人が来店していれば（N_t = 3）、X_t = U_1 + U_2 + U_3 である。足し合わせる個数 N_t そのものがランダムである点が、この確率過程の特徴である。

以下では、U_k の平均を \mu = E[U_k]、分散を \sigma^2 = V[U_k] と書く。パン屋の例では U_k は購入額、\lambda は1時間あたりの平均客数である。特に、すべての k で U_k = 1 なら X_t = N_t となり、ポアソン過程そのものに戻る。つまり複合ポアソン過程は、ポアソン過程の「出来事が1回起こるたびに1増える」を「出来事が1回起こるたびに U_k だけ増える」に置き換えたものである。

図1は、パン屋の例で1時間分のパスを描いたものである（購入額は平均1.5千円の指数分布に従う乱数とした）。上の N_t と下の X_t は同じ時刻に跳ね上がるが、跳ね上がる高さが違う。N_t は客が来るたびに1ずつ増えるのに対し、X_t はその客の購入額 U_k だけ増える。この例では1時間に8人が来店し、売上の合計は約11.7千円であった。

独立定常増分をもつ

複合ポアソン過程も、ポアソン過程と同じく独立定常増分をもつ。時刻 t から t + h までの増分 X_{t+h} - X_t は、その間に来店した客の購入額の合計である。その人数 N_{t+h} - N_t は過去と独立に Po(\lambda h) に従い、足し合わせる購入額も時刻 t より前の客のものとは別なので、次のことがわかる。

• 重ならない時間帯の売上は互いに独立である
• 長さ h の時間帯の売上 X_{t+h} - X_t は、時間帯の位置によらず X_h と同じ分布に従う

たとえば「10時台の売上」と「14時台の売上」は独立で、同じ分布に従う。この性質は、後の正規近似やパラメータ推定で使う。

平均と分散

売上の合計 X_t の平均と分散を求めよう。人数が n 人と決まっていれば、U_1 + \cdots + U_n の平均は n\mu、分散は（U_k が独立なので）n\sigma^2 である。ところが X_t では足し合わせる人数 N_t 自体がランダムなので、この公式をそのまま使うことはできない。そこで、次の2段階に分けて考える。

小さな例で考え方を確かめる

この2段階の計算を、人数の分布を単純にした例で実際にやってみる。開店直後の短い時間帯に来る客の数 N が0人、1人、2人のいずれかで、その確率がそれぞれ0.5、0.3、0.2だとする。1人あたりの購入額の平均は1500円で、人数とは無関係に決まるとする。

(1) の段階では、人数ごとに売上の平均を求める（表の右端の列）。(2) の段階では、それらを人数の確率で重みをつけて平均する。

一方、人数の平均は E[N] = 0 \times 0.5 + 1 \times 0.3 + 2 \times 0.2 = 0.7（人）なので、E[X] = 1500 \times 0.7 となっている。売上の平均は「1人あたりの平均額 × 平均人数」という自然な答えになった。

(1) で求めた「N = n のときの X の平均」を条件付き期待値といい、E[X \mid N = n] と書く。この例では E[X \mid N = n] = 1500n である。n に確率変数 N を入れたもの 1500N を E[X \mid N] と書くと、(2) の計算は E[X \mid N] の平均をとることにあたる。まとめると次のようになる。

これを繰り返し期待値の法則（全期待値の法則、law of total expectation）という。「場合分けして平均を求め、場合ごとの確率で平均し直す」という手順を式で書いたものである。

平均の導出

複合ポアソン過程でも、同じ2段階で計算する。

手順1（人数を固定する）　N_t = n のとき、X_t = U_1 + \cdots + U_n は n 個の和である。U_k は N と独立なので、人数がわかっても U_k の平均は \mu のままである。よって

手順2（人数について平均する）　n に N_t を入れると E[X_t \mid N_t] = \mu N_t である。これを N_t について平均すると、E[N_t] = \lambda t より

となる。パン屋の例では、1時間の売上の平均は E[X_1] = 8 \times 1 \times 1.5 = 12（千円）である。

分散の導出

分散も同じ2段階で求める。分散は V[X] = E[X^2] - (E[X])^2 で計算できるので、まず E[X_t^2] を求める。

手順1（人数を固定する）　N_t = n のとき、X_t = U_1 + \cdots + U_n は独立な n 個の和なので、平均は n\mu、分散は n\sigma^2 である。E[Y^2] = V[Y] + (E[Y])^2 の関係を使うと

手順2（人数について平均する）　n に N_t を入れて平均をとる。N_t \sim Po(\lambda t) より E[N_t] = V[N_t] = \lambda t なので、E[N_t^2] = V[N_t] + (E[N_t])^2 = \lambda t + (\lambda t)^2 である。よって

手順3（分散を計算する）　平均の2乗 (\lambda t \mu)^2 = \mu^2 (\lambda t)^2 を引くと、\mu^2 (\lambda t)^2 の項が打ち消し合って

最後の等号は \sigma^2 + \mu^2 = V[U_k] + (E[U_k])^2 = E[U_k^2] による。平均は「平均回数 × 1回あたりの平均」、分散は「平均回数 × 1回あたりの2乗の平均」と覚えておくとよい。

分散の2つの成分

V[X_t] = \lambda t \sigma^2 + \lambda t \mu^2 と2つに分けると、売上がばらつく2つの原因が見えてくる。

• \lambda t \sigma^2：購入額のばらつきによる部分。人数がちょうど平均どおりの \lambda t 人だったとしても、1人ごとの購入額が分散 \sigma^2 でばらつくので、合計の分散は \lambda t \sigma^2 になる。
• \lambda t \mu^2：人数のばらつきによる部分。全員がちょうど \mu ずつ買ったとしても、売上は \mu N_t となり、人数 N_t の分散 \lambda t によって分散 \mu^2 \lambda t が生じる。

複合ポアソン過程の分散は、この2つの和になっている。これは、一般に成り立つ全分散の公式（law of total variance）

の一例である。V[X \mid N = n] は人数を固定したときの分散（条件付き分散）で、ここでは E\bigl[V[X_t \mid N_t]\bigr] = E[\sigma^2 N_t] = \lambda t \sigma^2 が第1項、V\bigl[E[X_t \mid N_t]\bigr] = V[\mu N_t] = \lambda t \mu^2 が第2項にあたる。

パン屋の例では、購入額が指数分布に従うので \sigma = \mu = 1.5 であり、1時間の売上の分散は

となる。標準偏差は6千円で、2つの原因がちょうど半分ずつ分散に寄与している。

図2は、パン屋の8時間分のパスを3本描き、平均 E[X_t] = 12t と、平均 \pm 2 \times 標準偏差（標準偏差は \sqrt{36t} = 6\sqrt{t}）の線を重ねたものである。平均は時間に比例して増え、ばらつきの幅は \sqrt{t} に比例して広がる。分散が時間に比例して増えるのは、独立定常増分性により、1時間ごとの売上の分散36が時間とともに積み重なっていくためである。8時間の売上は平均96千円、標準偏差 6\sqrt{8} \approx 17.0 千円となる。

モーメント母関数

平均と分散だけでなく、X_t の分布そのものを調べるにはモーメント母関数が役に立つ。U_k のモーメント母関数を

とおく。X_t のモーメント母関数 M(s) = E\bigl[e^{sX_t}\bigr] も、平均と同じ2段階で求められる。

手順1（人数を固定する）　N_t = n のとき e^{sX_t} = e^{sU_1} e^{sU_2} \cdots e^{sU_n} であり、U_1, \dots, U_n は独立なので、期待値は積に分かれる。

手順2（人数について平均する）　P(N_t = n) = e^{-\lambda t} \dfrac{(\lambda t)^n}{n!} で重みをつけて足し合わせる。

2行目から3行目への変形では、指数関数の展開

を x = \lambda t \, \phi(s) として使った。

U_k = 1 のときは \phi(s) = e^s なので M(s) = \exp\{\lambda t (e^s - 1)\} となり、ポアソン分布 Po(\lambda t) のモーメント母関数に一致する。また手順2の計算は、ポアソン分布の確率母関数 G(z) = E\bigl[z^{N_t}\bigr] = e^{\lambda t (z - 1)} に z = \phi(s) を代入したものと見ることもできる。

モーメント母関数から平均と分散を確かめる

M(s) を微分すると、先ほど求めた平均と分散をもう一度確かめられる。\phi(0) = 1、\phi'(0) = E[U_k] = \mu、\phi''(0) = E[U_k^2]、M(0) = 1 に注意する。合成関数の微分により

となるので、s = 0 を代入すると

よって V[X_t] = M''(0) - M'(0)^2 = \lambda t \, E[U_k^2] となり、条件付き期待値で求めた結果と一致する。

間引き：ベルヌーイ変数を足し合わせる場合

複合ポアソン過程の重要な特別な場合として、U_k が0か1の値をとるベルヌーイ分布に従う場合がある。たとえば、Webサイトへの訪問が強度 \lambda のポアソン過程に従い、訪問者1人1人がほかの訪問者とは独立に確率 q で商品を購入するとしよう。k 番目の訪問者が購入すれば U_k = 1、購入しなければ U_k = 0 とすると、X_t = \sum_{k=1}^{N_t} U_k は時刻 t までの購入件数になる。このように、起こったイベントの中から一部だけを確率的に選び出すことを間引き（thinning）という。

このとき U_k のモーメント母関数は \phi(s) = (1 - q) \cdot e^{0} + q \cdot e^{s} = 1 - q + qe^s なので、X_t のモーメント母関数は

となる。これは Po(\lambda q t) のモーメント母関数そのものである。モーメント母関数と分布は1対1に対応するので、X_t \sim Po(\lambda q t) がわかる。さらに X は独立定常増分をもつので、X 自体が強度 \lambda q のポアソン過程になる。

平均と分散の公式からも確かめられる。U_k は0か1なので U_k^2 = U_k であり、\mu = E[U_k] = q、E[U_k^2] = q となる。よって E[X_t] = V[X_t] = \lambda q t で、平均と分散が等しいというポアソン分布の特徴と一致する。なお、データから q を推定するには、選ばれたイベントの強度の推定値を全体の強度の推定値 \hat{\lambda} で割ればよい。同じ時間だけ観測したなら、これは「購入件数 ÷ 訪問者数」に等しい。

なお、選ばれなかったイベント（購入しなかった訪問）も強度 \lambda (1 - q) のポアソン過程になり、しかも選ばれたイベントの過程とは互いに独立になることが知られている。1つのポアソン過程を確率的に振り分けると、互いに独立な2つのポアソン過程に分かれるのである。

計算例：保険金の支払い総額

導入で紹介した保険の例を計算してみよう。ある保険会社には、保険金の請求が強度 \lambda = 2（件/日）のポアソン過程に従って届く。1件あたりの請求額は平均 \mu = 10（万円）、標準偏差 \sigma = 5（万円）で、件数とは独立とする。30日間の支払い総額を X_{30}（万円）とする。

平均と標準偏差

30日間の請求件数の平均は \lambda t = 2 \times 30 = 60（件）なので

標準偏差は \sqrt{7500} \approx 86.6（万円）である。分散7500のうち、請求額のばらつきによる部分は \lambda t \sigma^2 = 1500、件数のばらつきによる部分は \lambda t \mu^2 = 6000 で、件数のばらつきのほうが4倍も大きい。

準備金が足りなくなる確率

保険会社が30日分の支払いのために750万円の準備金を用意したとする。支払い総額が準備金を超える確率 P(X_{30} > 750) を求めたい。

X_{30} は1日ごとの支払いの和 X_{30} = (X_1 - X_0) + (X_2 - X_1) + \cdots + (X_{30} - X_{29}) であり、独立定常増分性より、各日の支払いは互いに独立に同じ分布に従う。したがって中心極限定理により、X_{30} は近似的に正規分布 N(600, \; 7500) に従う。標準正規分布に従う確率変数を Z とすると

となる。

この近似の精度を確かめるため、請求額が平均10、分散25のガンマ分布 Ga(4, \; 2.5) に従う場合の正確な分布を図4に示す。正確な確率も、平均を求めたときと同じく件数で場合分けして計算できる。請求が n 件のとき、請求額の合計はガンマ分布の再生性より Ga(4n, \; 2.5) に従うので、その確率を件数の確率で重みをつけて足し合わせる。

和の計算は数値的に行った（Python実装を参照）。正規近似の0.042は、正確な値0.046をやや下回る。図4のとおり X_{30} の分布は右に少し裾が長く、左右対称な正規分布ではその分だけ右裾の確率が小さく見積もられるためである。それでも、おおよその大きさをつかむには正規近似で十分に役に立つ。

パラメータ推定

複合ポアソン過程のパラメータは、イベントの強度 \lambda と、大きさ U_k の分布（平均 \mu や分散 \sigma^2 など）である。どのようなデータが手に入るかで推定の方法が変わる。

(I) 1件ごとの時刻と大きさがわかる場合

レジの記録から、客1人1人の来店時刻と購入額がわかる場合である。来店時刻だけを見ればポアソン過程のデータであり、購入額 U_1, \dots, U_n は互いに独立に同じ分布に従う標本である。したがって、\lambda はポアソン過程のパラメータ推定と同じ方法で、\mu は購入額の標本平均で推定すればよい。

N と U_k は独立なので、尤度は「来店時刻の部分（\lambda だけを含む）」と「購入額の部分（U_k の分布のパラメータだけを含む）」の積に分かれる。そのため最尤推定でも、2つの部分を別々に最大化すればよい。1時間ごとの客数と売上が両方記録されている場合も同様で、\lambda は総客数を観測時間で割って、\mu は総売上を総客数で割って推定できる。

(II) 一定時間ごとの合計しかわからない場合

次に、1時間ごとの売上の合計 Y_1, \dots, Y_n しか記録されていない場合を考える。Y_j = X_{j\Delta} - X_{(j-1)\Delta}（\Delta = 1 時間）は、独立定常増分性より互いに独立に同じ分布に従い、平均と分散は

である。モーメント法では、これらを標本平均 m と標本分散 v（n で割ったもの）に等しいとおいて解く。ただし未知数は \lambda, \mu, \sigma^2 の3つで式は2本しかないので、U_k の分布に仮定をおく必要がある。ここでは購入額が指数分布に従うと仮定する。平均 \mu の指数分布の標準偏差は \mu なので \sigma^2 = \mu^2 となり

となる。2つの式の比をとると \dfrac{v}{m} = 2\mu なので、次の推定量が得られる。

計算例

パン屋のある日の1時間ごとの売上（千円）が、10時間分で次のようであったとする。

標本平均と標本分散は

である。よって

売上の合計しかわからなくても、平均と分散の関係から「1時間に約8人、1人あたり約1500円」と、客数と客単価を分けて推定できた。分散 \lambda \Delta (\sigma^2 + \mu^2) が、人数のばらつきと購入額のばらつきの両方を含んでいるおかげである。

練習問題

まとめ

Python実装

複合ポアソン過程を生成する関数と、平均・分散の公式、合計のデータからのモーメント法による推定を実装し、本文の例を確かめる。

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

def simulate_compound_poisson(lam, T, jump_sampler, rng):
    """強度 lam のポアソン過程の到着時刻と、各イベントの大きさ U_k を区間 [0, T] で生成する"""
    times = []
    t = rng.exponential(1 / lam)              # 到着間隔は平均 1/lam の指数分布
    while t <= T:
        times.append(t)
        t += rng.exponential(1 / lam)
    jumps = jump_sampler(len(times), rng)     # 各イベントの大きさ U_1, U_2, ...
    return np.array(times), jumps

def sample_totals(lam, t, jump_sampler, size, rng):
    """時刻 t での値 X_t = U_1 + ... + U_{N_t} を size 個生成する"""
    N = rng.poisson(lam * t, size)            # 各回の件数 N_t
    U = jump_sampler(N.sum(), rng)            # すべての回の U_k をまとめて生成
    owner = np.repeat(np.arange(size), N)     # 各 U_k が何回目のものか
    return np.bincount(owner, weights=U, minlength=size)

def mean_var(lam, t, mu, sigma2):
    """E[X_t] = λtμ,  V[X_t] = λt(σ² + μ²)"""
    return lam * t * mu, lam * t * (sigma2 + mu**2)

def estimate_exp_jumps(totals, dt):
    """間隔 dt ごとの合計から、U_k を指数分布と仮定してモーメント法で (λ, μ) を推定する"""
    m, v = np.mean(totals), np.var(totals)    # 標本平均と標本分散（n で割る）
    mu_hat = v / (2 * m)
    return m / (dt * mu_hat), mu_hat

rng = np.random.default_rng(3)

# ========== パン屋：λ = 8 人/時, 購入額は平均 1.5 千円の指数分布 ==========
print("【パン屋の売上（λ = 8 人/時, 購入額の平均 1.5 千円）】")
exp_jumps = lambda n, rng: rng.exponential(1.5, n)
times, jumps = simulate_compound_poisson(8, 1, exp_jumps, rng)
print(f"1時間のパスの例: 客数 {len(times)} 人, 売上 {jumps.sum():.2f} 千円")
print(f"  最初の3人の来店時刻（分）: {np.round(times[:3] * 60, 1)}, 購入額（千円）: {np.round(jumps[:3], 2)}")
E, V = mean_var(8, 1, 1.5, 1.5**2)
X1 = sample_totals(8, 1, exp_jumps, 100000, rng)
print(f"1時間の売上 理論値: 平均 {E:.2f}, 分散 {V:.2f}")
print(f"  シミュレーション（10万回）: 平均 {X1.mean():.2f}, 分散 {X1.var():.2f}")

# ========== 保険金：λ = 2 件/日, 請求額は Ga(4, 2.5)（平均 10, 分散 25） ==========
print("\n【保険金の支払い総額（30日間）】")
gamma_jumps = lambda n, rng: rng.gamma(4, 2.5, n)
E, V = mean_var(2, 30, 10, 25)
print(f"平均 {E:.0f}, 分散 {V:.0f}, 標準偏差 {np.sqrt(V):.2f}")
print(f"P(X_30 > 750) 正規近似: {norm.sf(750, E, np.sqrt(V)):.4f}")
# 件数 n で場合分けして足し合わせる（n 件の請求額の和は Ga(4n, 2.5)）
p_exact = sum(poisson.pmf(n, 60) * gamma.sf(750, a=4 * n, scale=2.5) for n in range(1, 200))
print(f"P(X_30 > 750) 正確な値: {p_exact:.4f}")
X30 = sample_totals(2, 30, gamma_jumps, 100000, rng)
print(f"P(X_30 > 750) シミュレーション（10万回）: {np.mean(X30 > 750):.4f}")
print(f"（誤り）件数を60件に固定した場合の正規近似: {norm.sf(750, E, np.sqrt(60 * 25)):.6f}")

# ========== 間引き：訪問 λ = 30 人/時, 購入確率 q = 0.1 ==========
print("\n【間引き（λ = 30 人/時, q = 0.1）】")
bern_jumps = lambda n, rng: (rng.random(n) < 0.1).astype(float)
P1 = sample_totals(30, 1, bern_jumps, 100000, rng)
print(f"1時間の購入件数: 平均 {P1.mean():.3f}, 分散 {P1.var():.3f}（理論値はどちらも λq = 3）")
print(f"購入が0件だった割合: {np.mean(P1 == 0):.4f}（理論値 e^-3 = {np.exp(-3):.4f}）")

# ========== パラメータ推定：1時間ごとの売上の合計だけから ==========
print("\n【パラメータ推定（購入額は指数分布と仮定）】")
sales = np.array([9, 17, 5, 12, 24, 10, 7, 20, 6, 10])     # 1時間ごとの売上（千円）
lam_hat, mu_hat = estimate_exp_jumps(sales, 1)
print(f"10時間分のデータ: λ の推定値 {lam_hat:.2f} 人/時, μ の推定値 {mu_hat:.2f} 千円")
long_sales = sample_totals(8, 1, exp_jumps, 1000, rng)        # 1000時間分のシミュレーション
lam_hat, mu_hat = estimate_exp_jumps(long_sales, 1)
print(f"1000時間分のシミュレーション: λ の推定値 {lam_hat:.2f}, μ の推定値 {mu_hat:.2f}（真の値 8, 1.5）")

1時間の売上のシミュレーションでは、平均と分散が理論値の12と36にほぼ一致している。保険の例では、正規近似の0.0416に対して、件数で場合分けした正確な値は0.0464、シミュレーションは0.0470となり、正規近似がやや小さめであることが確認できる。件数を固定した誤った計算では0.000054となり、危険を大きく過小評価してしまう。間引きの例では、1時間の購入件数の平均と分散がどちらも約3で、購入が0件だった割合も e^{-3} に近く、強度3のポアソン過程になっていることがわかる。最後のパラメータ推定では、1000時間分のシミュレーションデータから真の値に近い推定値が得られている。