電球が切れたら新しいものに交換する、というシンプルな状況を考えましょう。各電球の寿命は確率的に異なりますが、同じ種類の電球を使い続けるとします。最初の電球の寿命が $X_1$、2番目の寿命が $X_2$、3番目が $X_3$、…と続きます。このとき、「時刻 $t$ までに何回交換が起きたか」や「長時間の平均交換率はいくらか」を分析するのが更新過程(renewal process)の理論です。
更新過程は、ポアソン過程の自然な一般化です。ポアソン過程では到着間隔が指数分布に従いますが、更新過程では任意の非負分布に従います。この一般化により、はるかに幅広い現実のシステムをモデル化できます。ポアソン過程が持つ便利な無記憶性は失われますが、長時間の漸近的な振る舞いには驚くほど普遍的な法則が成り立ちます。
更新過程を理解すると、以下のような応用が可能です。
- 信頼性工学: 機器の交換間隔の分析と最適交換政策の設計に使われます。寿命分布がワイブル分布に従う部品の交換計画が典型例です
- 保険数理: 保険金請求の到着過程のモデル化に利用されます。請求額と請求間隔を組み合わせた複合更新過程が基本モデルです
- 在庫管理: 需要の発生パターンの分析に応用され、最適発注量の決定に活用されます
- 確率過程の一般理論: マルコフ連鎖の再帰性の解析、生存時間分析、待ち行列理論など、多くの分野の基礎ツールです
- 通信工学: パケット到着の統計的モデルとして、ネットワーク容量の設計に利用されます
本記事の内容
- 更新過程の定義と更新関数の性質
- 初等更新定理の証明と具体例
- 更新過程の中心極限定理
- ブラックウェルの更新定理と格子・非格子の区別
- キー更新定理と更新方程式の解法
- 更新報酬定理と応用
- 残余寿命・年齢・検査のパラドックスの完全な解析
- 遅延更新過程と定常更新過程
- Pythonによるシミュレーションと検証
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ポアソン分布 — ポアソン過程の基礎
- 指数分布 — 更新過程の特殊ケース
- 大数の法則と中心極限定理 — 漸近定理の基礎
更新過程の定義
基本的な構成
更新過程は、非負の独立同一分布(i.i.d.)の確率変数の列 $X_1, X_2, X_3, \dots$ から構成されます。$X_n$ は第 $n$ 回目の更新間隔(inter-renewal time)と呼ばれます。
電球の例に戻ると、$X_n$ は $n$ 番目の電球の寿命です。寿命は確率的に変動しますが、各電球は同じ種類なので統計的な性質(分布)は同じです。「独立」の仮定は、ある電球の寿命が他の電球の寿命に影響しないことを意味しています。
$X_n$ の共通分布関数を $F(x) = P(X_n \leq x)$、期待値を $\mu = E[X_n]$、分散を $\sigma^2 = \text{Var}(X_n)$ とします。技術的な仮定として $F(0) = P(X_n = 0) < 1$(更新間隔がほぼ確実に正)と $0 < \mu < \infty$(有限の平均寿命)を置きます。$F(0) < 1$ は「電球が即座に切れる確率が1未満」を意味し、これがないと有限時間に無限回の更新が起きうるという病理的な状況が生じます。
$n$ 回目の更新が起きる時刻(更新時刻、renewal epoch)は、更新間隔の部分和です。
$$ \begin{equation} S_n = X_1 + X_2 + \cdots + X_n, \quad S_0 = 0 \end{equation} $$
$S_n$ は $n$ 番目の電球が切れる時刻、あるいは $n$ 回目の交換が行われる時刻です。大数の法則により $S_n/n \to \mu$(a.s.)なので、$S_n \to \infty$(a.s.)であり、更新時刻は確率1で無限遠に発散します。
更新計数過程
時刻 $t$ までに起きた更新の回数を更新計数過程(renewal counting process)と呼びます。
$$ \begin{equation} N(t) = \max\{n \geq 0 : S_n \leq t\} = \sum_{n=1}^{\infty}\mathbf{1}(S_n \leq t) \end{equation} $$
$N(t)$ は右連続な階段関数で、各更新時刻 $S_n$ で1だけジャンプします。$N(t) = n$ は「時刻 $t$ までにちょうど $n$ 回の更新が起きた」($S_n \leq t < S_{n+1}$)を意味します。
$N(t)$ と $S_n$ には次の重要な双対関係があります。
$$ N(t) \geq n \iff S_n \leq t $$
この関係は、更新計数が更新時刻の「逆関数」のようなものであることを示しています。多くの証明でこの双対性が活用されます。
$S_n \to \infty$(a.s.)であることから $N(t) < \infty$(a.s.)が保証されます。つまり有限時間に無限回の更新が起きることはありません。さらに $N(t) \to \infty$($t \to \infty$)が成り立ちます。
ポアソン過程との関係
$X_n \sim \text{Exp}(\lambda)$(パラメータ $\lambda$ の指数分布)の場合、更新過程はパラメータ $\lambda$ のポアソン過程になります。このとき $N(t) \sim \text{Poisson}(\lambda t)$ です。
ポアソン過程が特別な理由は、指数分布の無記憶性(memoryless property)$P(X > s + t \mid X > s) = P(X > t)$ にあります。この性質により、「現在から次の更新までの待ち時間」の分布が、「前回の更新からどれだけ時間が経ったか」に依存しません。一般の更新過程では無記憶性が成り立たないため、過程の振る舞いがより複雑になります。
具体的に言うと、ポアソン過程では任意の時刻 $t$ から見て、次のイベントまでの待ち時間は常に $\text{Exp}(\lambda)$ に従います。しかし、例えばガンマ分布の更新間隔を持つ更新過程では、前回の更新から時間が経つにつれて「次の更新がもうすぐ来る」と予測できるようになり、残余寿命の分布が変化します。
この違いにもかかわらず、長時間の漸近的な性質においては、更新間隔の分布の詳細によらない普遍的な結果が成り立ちます。これが更新理論の魅力です。
更新関数
更新関数(renewal function)$m(t)$ は、時刻 $t$ までの更新回数の期待値です。
$$ \begin{equation} m(t) = E[N(t)] = \sum_{n=1}^{\infty}P(S_n \leq t) = \sum_{n=1}^{\infty}F^{*n}(t) \end{equation} $$
ここで $F^{*n}$ は分布関数 $F$ の $n$ 重畳み込みです。$F^{*1} = F$, $F^{*2}(t) = \int_0^t F(t-s) \, dF(s)$, … と再帰的に定義されます。$F^{*n}(t) = P(S_n \leq t)$ は $n$ 回目の更新が時刻 $t$ 以前に起きる確率です。
更新関数は更新過程の最も基本的な量であり、多くの理論的結果がこの関数を通じて表現されます。しかし、具体的な分布に対して $m(t)$ を閉じた形で求めることは通常困難です。ポアソン過程($m(t) = \lambda t$)は例外的に簡単なケースです。
$m(t)$ が任意の有限な $t$ に対して有限であること、すなわち $m(t) < \infty$($\forall t < \infty$)は、$F(0) < 1$ の仮定から証明できます。
ワルドの等式と更新関数の初等的評価
更新過程に関連する重要な結果として、ワルドの等式(Wald’s equation)があります。$N(t) + 1$ は停止時刻(正確には $N(t)+1$ は更新間隔列に対する停止時刻)であり、ワルドの等式から
$$ E[S_{N(t)+1}] = \mu \cdot E[N(t) + 1] = \mu \cdot (m(t) + 1) $$
$S_{N(t)+1} > t$(次の更新は $t$ より後に起きる)なので、$E[S_{N(t)+1}] > t$ です。よって、
$$ \mu \cdot (m(t) + 1) > t \implies m(t) > \frac{t}{\mu} – 1 $$
一方、$S_{N(t)} \leq t$ から $E[S_{N(t)}] \leq t$、すなわち $\mu \cdot m(t) \leq t$ より $m(t) \leq t/\mu$… と行きたいところですが、この不等式は $E[S_{N(t)}] = \mu \cdot E[N(t)]$ と直接は言えません($N(t)$ が確率変数なのでワルドの等式の適用条件に注意が必要です)。
実際、$S_{N(t)+1} = t + \gamma(t)$($\gamma(t)$ は残余寿命)なので、
$$ m(t) = \frac{t + E[\gamma(t)]}{\mu} – 1 $$
これは $m(t)$ と残余寿命の関係を示す基本的な等式です。
更新過程の基本概念が定義できたところで、最も重要な漸近結果である初等更新定理を証明しましょう。
初等更新定理
定理の主張
初等更新定理(elementary renewal theorem)は、更新過程の最も基本的な漸近結果です。2つの形があります。
強法則形(概収束):
$$ \begin{equation} \lim_{t \to \infty}\frac{N(t)}{t} = \frac{1}{\mu} \quad \text{a.s.} \end{equation} $$
更新関数形:
$$ \begin{equation} \lim_{t \to \infty}\frac{m(t)}{t} = \frac{1}{\mu} \end{equation} $$
つまり、長時間の平均では、単位時間あたりの更新回数は $1/\mu$(更新間隔の期待値の逆数)に収束します。
この結果は非常に直感的です。更新間隔の平均が $\mu$ なら、単位時間あたり平均 $1/\mu$ 回の更新が起きるはずです。例えば電球の平均寿命が1000時間なら、長期的には1時間あたり $1/1000 = 0.001$ 回の交換が期待されます。
証明(強法則形)
証明は大数の法則を巧みに利用します。$N(t)$ 回の更新が時刻 $t$ までに起きたとき、更新時刻について次の不等式が成り立ちます。
$$ S_{N(t)} \leq t < S_{N(t)+1} $$
左の不等式は「$N(t)$ 回目の更新は $t$ 以前に起きた」、右の不等式は「$N(t)+1$ 回目の更新は $t$ より後に起きる」ことを意味しています。
両辺を $N(t)$($N(t) > 0$ のとき)で割ると、
$$ \frac{S_{N(t)}}{N(t)} \leq \frac{t}{N(t)} < \frac{S_{N(t)+1}}{N(t)} $$
$N(t) \to \infty$(a.s.)なので、大数の法則 $S_n/n \to \mu$(a.s.)を $n = N(t)$ に適用すると($N(t) \to \infty$ かつ $N(t)$ は非減少なので適用可能)、
左辺: $\frac{S_{N(t)}}{N(t)} \to \mu \quad \text{(a.s.)}$
右辺: $\frac{S_{N(t)+1}}{N(t)} = \frac{S_{N(t)+1}}{N(t)+1} \cdot \frac{N(t)+1}{N(t)} \to \mu \cdot 1 = \mu \quad \text{(a.s.)}$
はさみうちの原理より $\frac{t}{N(t)} \to \mu$(a.s.)、つまり $\frac{N(t)}{t} \to \frac{1}{\mu}$(a.s.)が得られます。$\blacksquare$
この証明は大数の法則を使っており、更新間隔分布の具体的な形に全く依存しない「普遍的」な結果であることが明確です。
具体例: 更新率の計算
いくつかの具体的な分布に対して更新率 $1/\mu$ を計算しましょう。
| 更新間隔分布 | $\mu$ | $\sigma^2$ | 更新率 $1/\mu$ |
|---|---|---|---|
| $\text{Exp}(\lambda)$ | $1/\lambda$ | $1/\lambda^2$ | $\lambda$ |
| $\text{Uniform}(0, 2a)$ | $a$ | $a^2/3$ | $1/a$ |
| $\text{Gamma}(\alpha, \beta)$ | $\alpha/\beta$ | $\alpha/\beta^2$ | $\beta/\alpha$ |
| $\text{Weibull}(k, \lambda)$ | $\lambda\Gamma(1 + 1/k)$ | $\lambda^2[\Gamma(1+2/k) – \Gamma^2(1+1/k)]$ | $1/[\lambda\Gamma(1+1/k)]$ |
| $\text{LogNormal}(\mu_L, \sigma_L^2)$ | $e^{\mu_L + \sigma_L^2/2}$ | $(e^{\sigma_L^2}-1)e^{2\mu_L+\sigma_L^2}$ | $e^{-\mu_L – \sigma_L^2/2}$ |
すべての場合で、長時間の平均更新率は $1/\mu$ です。分布の形状(分散、歪度、尖度など)は漸近的な更新率には影響しません。しかし、有限時間での $m(t)$ の振る舞いや、更新回数の変動(標準偏差)には分散が影響します。
更新関数形の証明
更新関数形 $m(t)/t \to 1/\mu$ の証明には、概収束から期待値の収束を導く必要があり、技術的に少し難しくなります。鍵となるのは $N(t)/t$ の一様可積分性を示すことです。
$N(t) + 1 \leq \sum_{n=1}^{\infty} \mathbf{1}(X_n \leq t) + 1$ のような評価を用いて、$E[N(t)^2] < \infty$ を示し、一様可積分性を経由して結論を得ます。詳細は技術的なので省略しますが、結果は直感に合う自然なものです。
更新過程の中心極限定理
定理の主張と導出
更新回数 $N(t)$ は、初等更新定理で $N(t) \approx t/\mu$ と近似されますが、この近似のまわりの変動はどの程度でしょうか。更新過程の中心極限定理がこの問いに答えます。
$$ \begin{equation} \frac{N(t) – t/\mu}{\sigma\sqrt{t}/\mu^{3/2}} \xrightarrow{d} N(0, 1) \quad (t \to \infty) \end{equation} $$
ここで $\sigma^2 = \text{Var}(X_n)$ は更新間隔の分散です。
この結果は、$N(t)$ の標準偏差が $\sigma\sqrt{t}/\mu^{3/2}$ のオーダーであることを示しています。変動は $\sqrt{t}$ に比例して増加しますが、$N(t)$ 自体は $t/\mu$ に比例して増加するため、相対的な変動(変動係数)は $\sigma/(\mu\sqrt{t/\mu}) = \sigma/(t\mu)^{1/2} \cdot \mu^{-1/2}$ で $t$ とともに0に収束します。
証明のスケッチ
証明の方針は、$S_n$ に対する通常の中心極限定理を $N(t)$ に「逆変換」することです。
$S_n$ に対する中心極限定理: $(S_n – n\mu)/(\sigma\sqrt{n}) \xrightarrow{d} N(0, 1)$
$N(t) \geq n \iff S_n \leq t$ の双対性を使って、$P(N(t) \leq n)$ を $P(S_n \geq t)$ に変換します。$n = t/\mu + z\sigma\sqrt{t}/\mu^{3/2}$ と置くと、
$$ P\left(\frac{N(t) – t/\mu}{\sigma\sqrt{t}/\mu^{3/2}} \leq z\right) = P(N(t) \leq n) = P(S_n \geq t) \to 1 – \Phi(-z) = \Phi(z) $$
厳密には $n$ が整数である必要があるため、離散化の処理が必要ですが、本質的なアイデアは上の通りです。
信頼区間と応用
中心極限定理を使うと、$N(t)$ の近似的な信頼区間を構成できます。95%信頼区間は
$$ N(t) \in \left[\frac{t}{\mu} – 1.96\frac{\sigma\sqrt{t}}{\mu^{3/2}}, \quad \frac{t}{\mu} + 1.96\frac{\sigma\sqrt{t}}{\mu^{3/2}}\right] $$
例えば、平均寿命 $\mu = 100$ 時間、標準偏差 $\sigma = 50$ 時間の電球を使い続けた場合、$t = 10000$ 時間後の交換回数は約 $100 \pm 1.96 \times 50\sqrt{10000}/100^{3/2} = 100 \pm 0.98$、つまり99回から101回程度と予測されます。
初等更新定理は大域的な平均を記述しますが、次に局所的な更新密度の収束を保証するより精密な結果を見ましょう。
ブラックウェルの更新定理
定理の主張
初等更新定理は $m(t)/t \to 1/\mu$ を主張しますが、更新関数の増分についてより精密な結果を与えるのがブラックウェルの更新定理(Blackwell’s renewal theorem)です。
更新間隔分布 $F$ が非格子(non-lattice, non-arithmetic)の場合、任意の $h > 0$ に対して、
$$ \begin{equation} \lim_{t \to \infty} [m(t + h) – m(t)] = \frac{h}{\mu} \end{equation} $$
つまり、十分時間が経った後は、長さ $h$ の時間区間 $[t, t+h]$ に含まれる更新の期待回数は $h/\mu$ に近づきます。
非格子分布と格子分布
ブラックウェルの定理を正しく理解するには、非格子分布と格子分布の区別が重要です。
分布 $F$ が格子分布(lattice distribution, arithmetic distribution)であるとは、ある $d > 0$ が存在して $P(X \in \{d, 2d, 3d, \dots\}) = 1$ を満たすことです。このような最大の $d$ をスパン(span)と呼びます。例えば、更新間隔が常に整数値をとる場合は格子分布であり、スパンは1(または1の約数)です。
格子分布でない分布を非格子分布と呼びます。連続分布はすべて非格子分布です。混合分布(連続成分と離散成分の混合)の場合は個別の判定が必要です。
初等更新定理との違い
初等更新定理と比較すると、ブラックウェルの定理の精密さが際立ちます。
初等更新定理は「$[0, t]$ 全体での平均密度が $1/\mu$ に収束する」という大域的な結果です。更新関数を$[0, t]$ の長さ $t$ で平均するため、初期の過渡的な振る舞いが薄まって消えます。
一方、ブラックウェルの定理は「$[t, t+h]$ という局所的な区間での密度も $h/\mu$ に収束する」という、より精密な局所的結果です。初期の過渡現象が完全に消え去り、任意の局所区間で更新が「均質化」することを保証しています。
格子分布の場合
更新間隔分布がスパン $d$ の格子分布の場合、ブラックウェルの定理は次の形を取ります。
$$ \lim_{n \to \infty} u_n = \frac{d}{\mu} $$
ここで $u_n = P(S_k = nd \text{ for some } k \geq 1)$ は、時刻 $nd$ で更新が起きる確率です。
この結果は、1次元対称ランダムウォークの再帰性の証明でも使われます。ランダムウォークの原点帰還確率 $u_{2n} \sim 1/\sqrt{\pi n}$ の和が発散することが再帰性の証明に使われましたが、これは更新定理の一般的な枠組みの特殊ケースです。
更新密度の収束
更新間隔分布が密度 $f$ を持つ場合、更新関数 $m(t)$ も密度(更新密度)$m'(t)$ を持ちます。ブラックウェルの定理の微分版として、
$$ \lim_{t \to \infty} m'(t) = \frac{1}{\mu} $$
が成り立ちます。更新密度の具体的な表現は、更新方程式
$$ m'(t) = f(t) + \int_0^t m'(t-s) f(s) \, ds $$
を満たします。この積分方程式はラプラス変換 $\hat{m}'(s) = \hat{f}(s)/(1 – \hat{f}(s))$ を用いて解けます。
ブラックウェルの定理は更新関数の局所的な振る舞いを記述しますが、これをさらに一般化したのがキー更新定理です。
キー更新定理
定理の主張
ブラックウェルの定理の一般化がキー更新定理(key renewal theorem)です。この定理は更新過程に関連する多くの量の漸近挙動を統一的に導くための強力なツールです。
$h(t)$ を直接リーマン積分可能(directly Riemann integrable, d.R.i.)な関数とすると、
$$ \begin{equation} \lim_{t \to \infty} \int_0^t h(t – s) \, dm(s) = \frac{1}{\mu}\int_0^{\infty} h(t) \, dt \end{equation} $$
ここで $m(s)$ は更新関数であり、積分は更新測度に関するスティルチェス積分です。
直接リーマン積分可能性は技術的な条件ですが、実用上は以下の条件を満たす関数は d.R.i. です。
- 非負で、直接リーマン可積分($\int_0^{\infty} h(t) \, dt < \infty$)
- $t \to \infty$ で $h(t) \to 0$
- 有界で単調減少(十分条件)
典型的な応用では $h(t)$ は確率($0 \leq h(t) \leq 1$)や密度関数のように自然に積分可能な関数であり、d.R.i. の条件は問題にならないことが多いです。
キー更新定理がブラックウェルの定理を含むこと
$h(t) = \mathbf{1}_{[0, h_0]}(t)$(区間 $[0, h_0]$ の指示関数)を代入すると、
$$ \int_0^t \mathbf{1}_{[0, h_0]}(t – s) \, dm(s) = m(t) – m(t – h_0) $$
キー更新定理から $m(t) – m(t – h_0) \to h_0/\mu$ が得られ、これはブラックウェルの定理そのものです。
更新方程式とその解法
更新過程に関連する量 $g(t)$ は、しばしば次の形の更新方程式(renewal equation)を満たします。
$$ g(t) = h(t) + \int_0^t g(t – s) \, dF(s) $$
ここで $h(t)$ は既知の関数、$F$ は更新間隔分布です。
更新方程式の直感的な意味は、first-step analysisから理解できます。最初の更新が時刻 $s$($0 < s \leq t$)で起きたとします。その確率は $dF(s)$ です。最初の更新後、過程はリスタート(再生)するので、時刻 $s$ から先の寄与は $g(t-s)$ です。一方、最初の更新が起きる前($[0, s]$)の寄与が $h(t)$ に含まれています。
更新方程式の一般解はスティルチェス畳み込みで表現されます。
$$ g(t) = h(t) + \int_0^t h(t – s) \, dm(s) $$
キー更新定理を適用すると、解の漸近挙動が得られます。
$$ \lim_{t \to \infty} g(t) = \frac{1}{\mu}\int_0^{\infty} h(t) \, dt $$
この一般的な結果は、更新過程に関連するあらゆる量の長時間挙動を統一的に求めるための強力な道具です。
ラプラス変換による解法
更新方程式はラプラス変換を使うと代数方程式に変換できます。$\hat{g}(s) = \int_0^{\infty} e^{-st} g(t) \, dt$ とすると、
$$ \hat{g}(s) = \hat{h}(s) + \hat{g}(s) \hat{f}(s) $$
$$ \hat{g}(s) = \frac{\hat{h}(s)}{1 – \hat{f}(s)} $$
具体的な分布に対して $\hat{f}(s)$ を代入し、逆ラプラス変換で $g(t)$ を求められます。
更新報酬定理と応用
更新報酬定理
更新方程式の重要な応用として、更新報酬定理(renewal-reward theorem)があります。
各更新区間 $n$ でランダムな「報酬」$R_n$ が得られるとします。$(X_n, R_n)$ は i.i.d. とします(同じ更新区間内で $X_n$ と $R_n$ は相関してもよい)。時刻 $t$ までの累積報酬を $C(t) = \sum_{n=1}^{N(t)} R_n$ とすると、
$$ \lim_{t \to \infty} \frac{C(t)}{t} = \frac{E[R]}{E[X]} = \frac{E[R]}{\mu} \quad \text{a.s.} $$
長期的な「報酬率」は、1回の更新あたりの期待報酬を更新間隔の期待値で割ったものです。
応用例1: 交互更新過程と可用度
交互更新過程(alternating renewal process)では、システムが「稼働中」(on)と「故障中」(off)を交互に繰り返します。稼働時間 $U_n$ と故障修理時間 $D_n$ が独立で、それぞれ分布 $F_U$, $F_D$ に従うとします。
更新報酬定理を適用しましょう。1サイクルの長さは $X_n = U_n + D_n$(期待値 $\mu = E[U] + E[D]$)であり、1サイクルの「稼働時間」を報酬 $R_n = U_n$(期待値 $E[R] = E[U]$)とすると、長期的な可用度(availability)は
$$ \begin{equation} A = \lim_{t \to \infty} \frac{\text{稼働時間}}{t} = \frac{E[U]}{E[U] + E[D]} \end{equation} $$
キー更新定理を直接適用しても同じ結果が得られます。時刻 $t$ でシステムが稼働中である確率 $P(\text{on at } t)$ は更新方程式を満たし、
$$ \lim_{t \to \infty} P(\text{稼働中 at } t) = \frac{E[U]}{E[U] + E[D]} = A $$
具体例: サーバーの平均稼働時間が200時間、平均故障修理時間が2時間の場合、可用度は $A = 200/(200+2) = 0.990$(99.0%)です。稼働時間分布や故障時間分布の形状に関わらず、期待値のみで可用度が決まるという驚くべき結果です。
応用例2: 年齢交換政策
部品の年齢が $T$ に達したら予防交換(コスト $c_p$)し、$T$ に達する前に故障したら故障交換(コスト $c_f > c_p$)するという年齢交換政策を考えます。
1サイクルの期待コストは $E[R(T)] = c_p \bar{F}(T) + c_f F(T)$($\bar{F} = 1 – F$)、1サイクルの期待長さは $E[X(T)] = \int_0^T \bar{F}(t) \, dt$ です。更新報酬定理から、単位時間あたりの期待コストは
$$ C(T) = \frac{c_p \bar{F}(T) + c_f F(T)}{\int_0^T \bar{F}(t) \, dt} $$
$C(T)$ を最小化する $T^*$ が最適交換年齢です。この最適化問題は信頼性工学の基本問題です。
残余寿命と年齢
3つの関連する量
時刻 $t$ において、更新過程に関連する3つの重要な量があります。
残余寿命(excess life, residual life) $\gamma(t) = S_{N(t)+1} – t$: 次の更新までの残り時間
年齢(age, current life) $\delta(t) = t – S_{N(t)}$: 前回の更新からの経過時間
全寿命(total life, spread) $\beta(t) = S_{N(t)+1} – S_{N(t)} = \gamma(t) + \delta(t)$: 現在の更新間隔の長さ
バスを待つ人の立場で考えると、$\delta(t)$ は「前のバスが発車してからの時間」、$\gamma(t)$ は「次のバスが来るまでの時間」、$\beta(t)$ は「前のバスから次のバスまでの間隔」です。これらの量は互いに $\beta(t) = \gamma(t) + \delta(t)$ で関係しています。
残余寿命の定常分布
キー更新定理を用いると、$\gamma(t)$ の定常分布が求まります。$P(\gamma(t) > x)$ は更新方程式を満たし、$h(t) = \bar{F}(t + x) = 1 – F(t + x)$ として、キー更新定理を適用すると、
$$ \lim_{t \to \infty} P(\gamma(t) > x) = \frac{1}{\mu}\int_0^{\infty} \bar{F}(u + x) \, du = \frac{1}{\mu}\int_x^{\infty} \bar{F}(u) \, du $$
定常残余寿命の密度関数は
$$ f_{\gamma_\infty}(x) = \frac{1 – F(x)}{\mu} = \frac{\bar{F}(x)}{\mu} $$
これは元の更新間隔分布を長さで重み付けしたサイズバイアス分布(size-biased distribution)です。なぜ元の分布と異なるのでしょうか。ランダムな時刻 $t$ は、長い更新間隔の中に「入りやすい」ため、長い間隔が過剰にサンプリングされるのです。長さ $x$ の間隔に入る確率は $x$ に比例するため、密度は $x \cdot f(x)$ に比例し、正規化すると $x f(x) / \mu$ になります。残余寿命の密度 $\bar{F}(x)/\mu$ はこのサイズバイアスの結果です。
年齢の定常分布
年齢 $\delta(t)$ の定常分布は残余寿命と同じです。
$$ \lim_{t \to \infty} P(\delta(t) > x) = \frac{1}{\mu}\int_x^{\infty} \bar{F}(u) \, du $$
これは対称性から直感的に理解できます。定常状態では、更新間隔の中のランダムな位置にいるので、「始まりからの距離」と「終わりまでの距離」は同じ分布を持ちます。
残余寿命の期待値と検査のパラドックス
定常状態での残余寿命の期待値を計算しましょう。
$$ E[\gamma_\infty] = \int_0^{\infty} P(\gamma_\infty > x) \, dx = \frac{1}{\mu}\int_0^{\infty}\int_x^{\infty} \bar{F}(u) \, du \, dx $$
順序を入れ替えて($x$ が $0$ から $u$、$u$ が $0$ から $\infty$)、
$$ E[\gamma_\infty] = \frac{1}{\mu}\int_0^{\infty} u \bar{F}(u) \, du = \frac{E[X^2]}{2\mu} = \frac{\mu^2 + \sigma^2}{2\mu} = \frac{\mu}{2} + \frac{\sigma^2}{2\mu} $$
ここで $E[X^2] = \mu^2 + \sigma^2$ を使いました。
この式から検査のパラドックス(inspection paradox, waiting time paradox, bus paradox)が導かれます。
$E[\gamma_\infty] = \mu/2 + \sigma^2/(2\mu) \geq \mu/2$
等号は $\sigma^2 = 0$(確定的な更新間隔)のときのみ成り立ちます。つまり、ランダムな時刻に「検査」したとき、次の更新までの待ち時間の期待値は、更新間隔の平均の半分以上です。
具体例: バスの到着間隔が平均10分のポアソン過程に従うとします。ランダムに到着した乗客の待ち時間の期待値はいくらでしょうか。指数分布では $\sigma^2 = \mu^2$ なので、
$$ E[\gamma_\infty] = \frac{\mu}{2} + \frac{\mu^2}{2\mu} = \frac{\mu}{2} + \frac{\mu}{2} = \mu = 10 \text{ 分} $$
「5分」ではなく「10分」です。これは指数分布の無記憶性の帰結です。いつバス停に到着しても、前のバスが出発してからの時間は関係なく、次のバスまでの待ち時間は常に平均10分の指数分布に従います。
各分布でのパラドックスの大きさを比較しましょう。
| 更新間隔分布 | $\mu$ | $\sigma^2$ | $E[\gamma_\infty]$ | $E[\gamma_\infty]/(\mu/2)$ |
|---|---|---|---|---|
| 確定的 | $\mu$ | 0 | $\mu/2$ | 1.0 |
| 一様 $U(0, 2\mu)$ | $\mu$ | $\mu^2/3$ | $2\mu/3$ | 1.33 |
| ガンマ$(2, 2/\mu)$ | $\mu$ | $\mu^2/2$ | $3\mu/4$ | 1.5 |
| 指数分布 | $\mu$ | $\mu^2$ | $\mu$ | 2.0 |
| パレート($\alpha=2.5$) | $\mu$ | $3\mu^2$ | $2\mu$ | 4.0 |
更新間隔の変動(分散)が大きいほど、検査のパラドックスの効果が強くなります。確定的な間隔(バスが正確に10分ごとに来る)では、平均待ち時間は直感通り5分ですが、指数分布(ランダムに来るバス)では10分に倍増します。パレート分布のような重い裾を持つ分布ではさらに大きくなります。
全寿命の定常分布
全寿命 $\beta(t) = \gamma(t) + \delta(t)$ の定常分布も求めましょう。
$$ \lim_{t \to \infty} P(\beta(t) > x) = \frac{1}{\mu}\int_0^x u \, dF(u) + \bar{F}(x) \cdot \frac{x}{\mu} $$
… を整理すると、密度は
$$ f_{\beta_\infty}(x) = \frac{x f(x)}{\mu} $$
これはサイズバイアス分布そのものです。長い更新間隔ほどランダムな時刻に「観測されやすい」という効果を直接反映しています。
$E[\beta_\infty] = E[X^2]/\mu = \mu + \sigma^2/\mu \geq \mu$ であり、全寿命の期待値は更新間隔の平均以上です。
遅延更新過程と定常更新過程
遅延更新過程
遅延更新過程(delayed renewal process)では、最初の更新間隔 $X_1$ の分布が $F_1$($F$ と異なりうる)であり、$X_2, X_3, \dots$ は通常通り i.i.d.(分布 $F$)です。
漸近的性質(初等更新定理、ブラックウェルの定理、キー更新定理)は通常の更新過程と同じですが、有限時間での $m(t)$ の振る舞いは $F_1$ によって異なります。遅延更新過程の更新関数は
$$ m_d(t) = F_1(t) + \int_0^t m(t – s) \, dF_1(s) $$
を満たします。ここで $m(t)$ は通常の(遅延なし)更新関数です。
定常更新過程
特に重要な遅延更新過程が、$X_1$ の分布を
$$ F_1(x) = F_e(x) = \frac{1}{\mu}\int_0^x [1 – F(u)] \, du $$
と選んだ場合の定常更新過程(equilibrium renewal process, stationary renewal process)です。
$F_e$ は平衡分布(equilibrium distribution)と呼ばれ、残余寿命の定常分布と一致します。
定常更新過程は、あたかも無限の過去から動いていた更新過程の「定常状態」を時刻0から観察したものに相当します。定常更新過程では更新関数が
$$ m_d(t) = \frac{t}{\mu} $$
正確に成り立ちます(漸近的にではなく全時刻で)。
定常更新過程のもう一つの重要な性質は、残余寿命 $\gamma(t)$ と年齢 $\delta(t)$ の分布が時刻 $t$ に依存しないことです。すべての $t \geq 0$ で、$\gamma(t)$ と $\delta(t)$ は定常残余寿命分布に従います。
Pythonによるシミュレーション
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
def simulate_renewal(dist_func, T_max, seed=42):
"""更新過程のシミュレーション"""
np.random.seed(seed)
times = [0.0]
t = 0.0
while t < T_max:
x = dist_func()
t += x
if t <= T_max:
times.append(t)
return np.array(times)
T_max = 1000
exp_times = simulate_renewal(lambda: np.random.exponential(1.0), T_max)
gamma_times = simulate_renewal(lambda: np.random.gamma(3, 1/3), T_max)
unif_times = simulate_renewal(lambda: np.random.uniform(0.5, 1.5), T_max)
fig, axes = plt.subplots(2, 2, figsize=(13, 10))
# N(t)/t の収束
ax = axes[0, 0]
for times, label in [(exp_times, "Exponential"),
(gamma_times, "Gamma(3, 1/3)"),
(unif_times, "Uniform(0.5, 1.5)")]:
t_grid = np.linspace(10, T_max, 200)
Nt = [np.sum(times <= t) / t for t in t_grid]
ax.plot(t_grid, Nt, linewidth=1.5, label=label, alpha=0.8)
ax.axhline(1.0, color="black", linestyle="--", linewidth=1, label="1/mu=1")
ax.set_xlabel("Time t"); ax.set_ylabel("N(t)/t")
ax.set_title("Elementary renewal theorem"); ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# サンプルパス
ax = axes[0, 1]
for times, label, color in [(exp_times, "Exp", "blue"),
(gamma_times, "Gamma", "red"),
(unif_times, "Uniform", "green")]:
mask = times <= 50
ax.step(times[mask], np.arange(np.sum(mask)), where="post",
linewidth=1.5, label=label, color=color)
ax.set_xlabel("Time t"); ax.set_ylabel("N(t)")
ax.set_title("Renewal counting process"); ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# 更新間隔ヒストグラム
ax = axes[1, 0]
for times, label, color in [(exp_times, "Exp", "blue"),
(gamma_times, "Gamma", "red"),
(unif_times, "Uniform", "green")]:
ax.hist(np.diff(times), bins=40, density=True, alpha=0.3,
color=color, edgecolor="white", label=label)
ax.set_xlabel("Inter-renewal time"); ax.set_ylabel("Density")
ax.set_title("Inter-renewal distributions"); ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# CLTの検証
ax = axes[1, 1]
n_sim, t_eval = 5000, 500
N_vals = np.zeros(n_sim)
for i in range(n_sim):
np.random.seed(i)
t, count = 0.0, 0
while t < t_eval:
t += np.random.exponential(1.0)
if t <= t_eval: count += 1
N_vals[i] = count
z = (N_vals - t_eval) / np.sqrt(t_eval)
ax.hist(z, bins=50, density=True, alpha=0.6, color="steelblue",
edgecolor="white", label="Simulated")
x_n = np.linspace(-4, 4, 200)
ax.plot(x_n, stats.norm.pdf(x_n), "r-", linewidth=2, label="N(0,1)")
ax.set_xlabel("Standardized N(t)"); ax.set_ylabel("Density")
ax.set_title(f"CLT for renewal process (t={t_eval})"); ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このシミュレーション結果から、以下の点が確認できます。
-
初等更新定理の収束(左上): 3つの異なる更新間隔分布(指数、ガンマ、一様)いずれにおいても、$N(t)/t$ は $1/\mu = 1.0$ に収束しています。初期の過渡期を過ぎると、分布の形状に関わらず漸近的な更新率は同じです。一様分布は分散が小さいため早く収束し、指数分布は分散が大きいため変動が大きいことも確認できます
-
サンプルパスの違い(右上): 指数分布(ポアソン過程)のサンプルパスはジャンプ間隔がばらついていますが、一様分布のサンプルパスはほぼ等間隔に近くなっています。ガンマ分布はその中間です。これは更新間隔の変動係数 $\sigma/\mu$ の違いを反映しています
-
更新間隔の分布(左下): 3つの分布の形状は大きく異なりますが、いずれも平均1で共通しており、初等更新定理が漸近的な更新率の普遍性を保証する根拠です
-
中心極限定理(右下): 標準化された $N(t)$ のヒストグラムが標準正規分布によく適合しており、更新過程の中心極限定理が数値的に確認されています
検査のパラドックスのシミュレーション
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
n_sim = 50000
T_obs = 100.0
residual_exp, residual_unif = [], []
for _ in range(n_sim):
t = 0.0
while True:
x = np.random.exponential(1.0)
if t + x > T_obs:
residual_exp.append(t + x - T_obs)
break
t += x
for _ in range(n_sim):
t = 0.0
while True:
x = np.random.uniform(0.5, 1.5)
if t + x > T_obs:
residual_unif.append(t + x - T_obs)
break
t += x
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
ax = axes[0]
ax.hist(residual_exp, bins=60, density=True, alpha=0.7,
color='steelblue', edgecolor='white', label='Simulation')
x_r = np.linspace(0, 5, 200)
ax.plot(x_r, np.exp(-x_r), 'r-', linewidth=2, label='Theory: Exp(1)')
ax.set_xlabel('Residual life'); ax.set_ylabel('Density')
ax.set_title(f'Exp renewals: E[gamma]={np.mean(residual_exp):.3f} (theory: 1.0)')
ax.legend(); ax.grid(True, alpha=0.3)
ax = axes[1]
ax.hist(residual_unif, bins=60, density=True, alpha=0.7,
color='steelblue', edgecolor='white', label='Simulation')
mu_u, sig2_u = 1.0, 1/12
E_th = mu_u/2 + sig2_u / (2*mu_u)
ax.axvline(np.mean(residual_unif), color='red', linestyle='--',
linewidth=2, label=f'E[gamma]={np.mean(residual_unif):.3f}')
ax.axvline(E_th, color='green', linestyle=':', linewidth=2,
label=f'Theory: {E_th:.3f}')
ax.set_xlabel('Residual life'); ax.set_ylabel('Density')
ax.set_title('Uniform renewals: residual life')
ax.legend(); ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のパネルでは、指数分布の更新過程における残余寿命が再び指数分布に従うことが確認できます。これは指数分布の無記憶性の現れであり、$E[\gamma_\infty] = \mu = 1$ です。「待ち時間の平均が間隔の平均と等しい」という検査のパラドックスがシミュレーションで裏付けられています。
右のパネルでは、一様分布 $U(0.5, 1.5)$ の更新過程での残余寿命を示しています。理論値 $E[\gamma_\infty] = \mu/2 + \sigma^2/(2\mu) = 0.5 + (1/12)/(2 \times 1) \approx 0.542$ とシミュレーション結果がよく一致しています。一様分布では分散が小さいため、検査のパラドックスの効果は指数分布ほど顕著ではありません。
まとめ
本記事では、更新過程の理論を初等更新定理からPythonでの検証まで解説しました。
- 更新過程はi.i.d.の非負確率変数列から構成され、ポアソン過程の一般化です。更新計数 $N(t)$ と更新時刻 $S_n$ の双対関係 $N(t) \geq n \iff S_n \leq t$ が理論の要です
- 更新関数 $m(t) = E[N(t)] = \sum F^{*n}(t)$ は更新過程の基本的な量であり、ラプラス変換 $\hat{m}'(s) = \hat{f}(s)/(1-\hat{f}(s))$ で特徴づけられます
- 初等更新定理は $N(t)/t \to 1/\mu$(a.s.)を主張し、大数の法則から導かれます。更新率は更新間隔の平均のみで決まる普遍的な結果です
- 更新過程のCLTは $N(t)$ の変動が $\sigma\sqrt{t}/\mu^{3/2}$ のオーダーであることを示します
- ブラックウェルの定理は局所的な更新密度の収束 $m(t+h) – m(t) \to h/\mu$ を保証し、初等更新定理より精密な結果です
- キー更新定理は更新過程の漸近解析の万能ツールであり、更新方程式の解の漸近挙動を $\frac{1}{\mu}\int_0^{\infty} h(t) \, dt$ で与えます
- 更新報酬定理は長期平均報酬率 $E[R]/\mu$ を与え、可用度計算や最適交換政策の設計に応用されます
- 検査のパラドックス: 残余寿命の期待値は $\mu/2 + \sigma^2/(2\mu) \geq \mu/2$ であり、更新間隔の変動が大きいほどパラドックスの効果が強くなります
- 定常更新過程: 平衡分布 $F_e$ で遅延させると、$m_d(t) = t/\mu$ が全時刻で正確に成り立ちます
次のステップとして、以下の記事も参考にしてください。
- 連続時間マルコフ連鎖 — CTMCの理論
- 出生死亡過程 — CTMCの重要な応用
- ランダムウォークの理論 — 離散的な確率的遍歴