宮沢政清著「対話・確率過程入門」は、「先生」と生徒たちふたりとの対話形式で行う、確率過程についての解説書です。
まだ全部読み終えてはいませんが、この初版本では、第２章20ページの下から３つめと２つめの等式が誤植になっていて、その部分についてノートしたものです。
この最後には、この誤植ページの主題である「指数分布」の性質をJavaScriptでプログラム実行することによって確認する話を載せておきました。
20ページの誤植部分についてのノート
まず第2章の内容は、バス待ちなどの待ち時間を題材にしたもので、(バスの)到着時刻を確率変数$T$として扱う話になっています。このとき、(自分がバス停に来た)時刻$x$から見て$T$で到着した場合の到着待ち時間は$T-x$であり、この到着待ち時間$T-x$の期待値について考えていく話です。
まず「先生」よって、この待ち時間の条件付き期待値は、その到着時間の密度関数を$f(x)$とおいて、
\mathbb{E}(T - x | T > x) = \frac{1}{\mathbb{P}(T > x)}\int_{x}^{\infty} (u-x)f(u)du
と示されています。個別の到着時刻$u$に対して、待ち時間の量が$u-x$、$u$の確率が$f(u)$であり、待ち時間の期待値としてこれらを掛けて総和(積分)しているのでしょう。
一方、条件として、到着時刻$T$は$x$以降の時刻のなかだけに限定($T > x$)されます。この条件となる、全到着時刻全域に対して、到着時刻$T$が$x$以降の時刻である場合の割合が、$\mathbb{P}(T > x)$です。
この時点では、分布や密度に具体的なものは与えられていません。
20ページでは、この待ち時間期待値が時刻$x$によらずに一定の値$\mathbb{E}(T)$になるとした場合の、到着時刻$T$の分布関数$F(x)$とはどのようなものか、について考える部分です。
16ページの$F(x) = \mathbb{P}(X \leq x)$より、$\mathbb{P}(T > x) = 1 - F(x)$となるのでしょう。
そして、期待値が$x$に依存しない条件なし期待値$\mathbb{E}(T)$の場合に対して考えることから、
\frac{1}{1-F(x)}\int_{x}^{\infty} (u-x)f(u)du = \mathbb{E}(T)
とする等式を設定し、これをもとに下から３つ目の等式が作られているでしょう。
そして、書籍での２０ページ下から３つ目の式は、
\int_{x}^{\infty} (u-x)f(u)dx = \mathbb{E}(T)(1-F(x))
となっています。これは単に
\int_{x}^{\infty} (u-x)f(u)du = \mathbb{E}(T)(1-F(x))
の誤植でしょう。
そして20ページでは、続けざまにこの両辺を微分したものとして、下から２つ目の式、
1 - F(x) = -\mathbb{E}(T)f(x)
が、いきなり「先生」から出されます。
しかし、これは実際に計算したものでは
-1 + F(x) = -\mathbb{E}(T)f(x)
となるもので、書籍中のとは符号が違っています。
そしてその次の一番下の等式では、
\frac{\frac{d}{dx}(1-F(x))}{1-F(x)} = -\frac{1}{\mathbb{E}(T)}
となっています。ここでこの等式をもとに、下から二番目の式の形にもどしても、
-(1 - F(x)) = \mathbb{E}(T)\frac{d}{dx}(1-F(x))
となるので、この面からも書籍での下から二番目の式が違っていることは確認できるでしょう。
また17ページの
f(x) = \frac{d}{dx}F(x)
から、
\frac{d}{dx}(1-F(x)) = -f(x)
となることで、この結果を下から二番目では右辺に適用されていますが、そこから導出する一番下の式で微分型に戻しているので、これを適用する必要もないのではと思いました。
そもそも2つ目の等式の左辺の微分は、書籍でのように「先生」が直接だして、それを生徒たちもすぐ受け入れられる単純さでもないと思いました。実際、読むことで理解させるレベルにするなら、
\begin{align}
\frac{d}{dx}\int_{x}^{\infty} (u-x)f(u)du &= \lim_{t \to \infty} \frac{d}{dx}\int_{x}^{t}(u-x)f(u)du \\
&= \lim_{t \to \infty}\Big(\frac{d}{dx}\int_{x}^{t}uf(u)du -\frac{d}{dx}\big(x\int_{x}^{t}f(u)du\big)\Big)\\
&= \lim_{t \to \infty}\Big(-xf(x) - \big(\int_{x}^{t}f(u)du + x(-f(x))\big)\Big)\\
&= \lim_{t \to \infty}[-F(u)]_{x}^{t}\\
&= \lim_{t \to \infty}(-F(t)+F(x))
\end{align}
くらいのステップは必要だと思います(2行目から3行目では、$\int_{a}^{b}g(t)dt = -\int_{b}^{a}g(t)dt$と、$\frac{d}{dx}\int_{a}^{x}g(t)dt = g(x)$を組み合わせて使っています)。この結果に16ページにある、
\lim_{x \to \infty}F(x) = 1
を使うことで、先の左辺の$-1+F(x)$が出くることになります。
話を途中で止めるのもなんなので、以下は、この20ページからの話の結末について紹介をしておきます。
この20ページの一番下の等式に対して、$\lambda = \frac{1}{\mathbb{E}(x)}$とおいた微分方程式$\frac{d}{dx}(1-F(x)) = -\lambda(1-F(x))$をとくと、$1-F(x) = Ce^{-\lambda x}$となります。
この結果へ、分布関数の条件$F(0)=0$を想定すると、
F(x) = 1 - e^{-\lambda x}
という分布関数が得られ、この分布が「指数分布」とよばれるものである、といった結末になっています。
指数分布の性質をプログラムで確認する
最後に、プログラミングの話題となるよう、プログラム実行によってこの時刻にかかわらず待ち時間期待値が同じである性質を確認してみます。ここではJavaScriptを使います。
通常の言語環境にビルトインされている乱数機能は一様分布です。
一様分布を用いてある分布関数に従う乱数値を作るには、その分布の逆関数へ0から1の範囲の一様乱数の値を渡すことで実現できます。分布関数はマイナス無限大から引数値までの範囲への累積確率で、0から1の範囲で結果を返すからです。
指数分布は$F(x) = 1-e^{-\lambda x}$なので、この逆関数は$F^{-1}(y) = \frac{-\log(1-y)}{\lambda}$となります。
よって、指数分布に従うランダムな到着時間を返す関数は、以下のようになります:
function invExpoDist(y, l) {
// y = 1-exp(-l*x) => x = -log(1-y)/l
return -Math.log(1 - y) / l;
}
function randArrive(ev) {
return invExpoDist(Math.random(), 1 / ev);
}
ここではランダムに到着時刻を返す関数randArrive()のパラメータとしては、$\lambda$ではなく、その逆数である待ち時間の期待値で与えるようにしてあります。
続いて、指数分布の到着時刻リストの内容を一覧するために、数値配列を対数ヒストグラムでコンソールにプロットするコードも用意しておきます。
function plotLogHist(vs, start, end, count, unit=2) {
const span = (end - start) / count;
const hist = [...Array(count + 1).keys()].map(
i => vs.filter(
v => start + span * i <= v && v < start + span * (i + 1)).length);
hist.forEach((len, i) => {
const bar = "*".repeat((Math.log(len) / Math.log(unit)) >>> 0).padEnd(40);
const g = `${start + i * span}`.padStart(6);
console.log(`${g} : ${bar} ... ${len}`);
});
}
指数分布に従う10000個の到着時間に対し、待ち時間の期待値expectValueとして、3, 30, 300を用いて、それぞれに15個の時刻xを起点とした、それぞれのxでの待ち時間の平均を計算して表示するコードは、以下のようになります。
{
const total = 10000, xcount = 15;
const evs = [3, 30, 300];
for (const expectValue of evs) {
const arrives = [...Array(total).keys()].map(_ => randArrive(expectValue));
const xs = [...Array(xcount).keys()].map(i => expectValue / 3 * i);
console.log(`\n[expected value: ${expectValue}]`);
plotLogHist(arrives, 0, expectValue * 10, 15);
for (const x of xs) {
const ts = arrives.filter(t => t > x);
const ave = ts.reduce((r, t) => r + (t - x), 0) / ts.length;
const xText = `${x}`.padEnd(6), samples = `${ts.length}`.padEnd(8);
console.log(`x = ${xText} samples = ${samples} wait-average = ${ave}`);
}
}
}
nodejsのREPLやブラウザのJavaScriptコンソールなどを使って、これら3つのコードを上から順に実行した結果は、以下のようになりました:
[expected value: 3]
0 : ************ ... 4828
2 : *********** ... 2473
4 : ********** ... 1295
6 : ********* ... 699
8 : ******** ... 374
10 : ******* ... 145
12 : ****** ... 84
14 : ***** ... 51
16 : **** ... 19
18 : *** ... 9
20 : *** ... 9
22 : *** ... 11
24 : ... 0
26 : * ... 3
28 : ... 0
30 : ... 0
x = 0 samples = 10000 wait-average = 3.036644482859655
x = 1 samples = 7152 wait-average = 3.055554064891443
x = 2 samples = 5172 wait-average = 3.047325374304088
x = 3 samples = 3741 wait-average = 3.030530029617039
x = 4 samples = 2699 wait-average = 3.0248490986075702
x = 5 samples = 1979 wait-average = 2.9583329467617134
x = 6 samples = 1404 wait-average = 2.966193030503541
x = 7 samples = 1001 wait-average = 2.966703148686848
x = 8 samples = 705 wait-average = 3.0103394156622683
x = 9 samples = 491 wait-average = 3.1115940135853926
x = 10 samples = 331 wait-average = 3.396370331853496
x = 11 samples = 242 wait-average = 3.451983799175673
x = 12 samples = 186 wait-average = 3.3442608431385206
x = 13 samples = 134 wait-average = 3.4585332344476036
x = 14 samples = 102 wait-average = 3.3838442463261282
[expected value: 30]
0 : ************ ... 4894
20 : *********** ... 2464
40 : ********** ... 1287
60 : ********* ... 653
80 : ******** ... 347
100 : ******* ... 180
120 : ****** ... 81
140 : ***** ... 39
160 : **** ... 28
180 : *** ... 13
200 : ** ... 6
220 : * ... 3
240 : * ... 2
260 : ... 1
280 : ... 1
300 : ... 0
x = 0 samples = 10000 wait-average = 30.077136083470045
x = 10 samples = 7132 wait-average = 30.26234196627789
x = 20 samples = 5106 wait-average = 30.37995893886442
x = 30 samples = 3685 wait-average = 30.216580483483984
x = 40 samples = 2642 wait-average = 30.21005675612523
x = 50 samples = 1891 wait-average = 30.340614901969964
x = 60 samples = 1355 wait-average = 30.3852641050147
x = 70 samples = 980 wait-average = 30.148698123825405
x = 80 samples = 702 wait-average = 30.239340893828317
x = 90 samples = 509 wait-average = 29.818656236799463
x = 100 samples = 355 wait-average = 30.80635266448583
x = 110 samples = 255 wait-average = 30.905986151061466
x = 120 samples = 175 wait-average = 32.76071069914236
x = 130 samples = 125 wait-average = 33.94591643602512
x = 140 samples = 94 wait-average = 33.79990460531498
[expected value: 300]
0 : ************ ... 4884
200 : *********** ... 2463
400 : ********** ... 1299
600 : ********* ... 642
800 : ******** ... 359
1000 : ******* ... 178
1200 : ****** ... 72
1400 : ***** ... 43
1600 : ***** ... 35
1800 : *** ... 12
2000 : ** ... 7
2200 : * ... 3
2400 : * ... 2
2600 : ... 0
2800 : ... 0
3000 : ... 0
x = 0 samples = 10000 wait-average = 300.1291269539622
x = 100 samples = 7118 wait-average = 302.6665365459592
x = 200 samples = 5116 wait-average = 302.68527557176986
x = 300 samples = 3693 wait-average = 300.80284284558604
x = 400 samples = 2653 wait-average = 300.12532444942
x = 500 samples = 1879 wait-average = 304.3857249230324
x = 600 samples = 1354 wait-average = 304.78413881002336
x = 700 samples = 998 wait-average = 296.58360453816886
x = 800 samples = 712 wait-average = 296.3437734497861
x = 900 samples = 495 wait-average = 304.9161466951072
x = 1000 samples = 353 wait-average = 308.6132624021511
x = 1100 samples = 240 wait-average = 333.39322737733613
x = 1200 samples = 175 wait-average = 338.29614443641884
x = 1300 samples = 133 wait-average = 327.7802449175222
x = 1400 samples = 103 wait-average = 311.13745973387717
この対数ヒストグラムは、1目盛りで数が2倍になっているものです。このため指数分布からのサンプルは、均一区間で一定比率で減少していく分布であることが確認できるでしょう。ちなみに、半減する区間は$\log(2)/\lambda$、つまり期待値にlog(2)=0.693...を掛けた値です。
時刻が後ろに行くほど対象となるサンプル数が減ってしまうので理想から外れていきますが、サンプル数が十分あることを想像すれば、どの時刻でも平均値が期待値に近いものとなることがわかるでしょう。