「対話・確率過程入門」初版の第２章の誤植を検証し、ついでに指数分布の性質をプログラム実行によって確認してみた

宮沢政清著「[対話・確率過程入門](https://www.gensu.jp/product/%E5%AF%BE%E8%A9%B1%E3%83%BB%E7%A2%BA%E7%8E%87%E9%81%8E%E7%A8%8B%E5%85%A5%E9%96%80/)」は、「先生」と生徒たちふたりとの対話形式で行う、確率過程についての解説書です。

まだ全部読み終えてはいませんが、この初版本では、第２章20ページの下から３つめと２つめの等式が誤植になっていて、その部分についてノートしたものです。

この最後には、この誤植ページの主題である「指数分布」の性質をJavaScriptでプログラム実行することによって確認する話を載せておきました。

## 20ページの誤植部分についてのノート

まず第2章の内容は、バス待ちなどの待ち時間を題材にしたもので、(バスの)**到着時刻を確率変数$T$**として扱う話になっています。このとき、(自分がバス停に来た)時刻$x$から見て$T$で到着した場合の到着待ち時間は$T-x$であり、この**到着待ち時間$T-x$の期待値**について考えていく話です。

まず「先生」よって、この待ち時間の条件付き期待値は、その到着時間の密度関数を$f(x)$とおいて、

```math
\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)$の場合に対して考えることから、

```math
\frac{1}{1-F(x)}\int_{x}^{\infty} (u-x)f(u)du = \mathbb{E}(T)
```

とする等式を設定し、これをもとに下から３つ目の等式が作られているでしょう。

---

そして、書籍での２０ページ下から３つ目の式は、

```math
\int_{x}^{\infty} (u-x)f(u)dx = \mathbb{E}(T)(1-F(x))
```

となっています。これは単に

```math
\int_{x}^{\infty} (u-x)f(u)du = \mathbb{E}(T)(1-F(x))
```

の誤植でしょう。

---

そして20ページでは、続けざまにこの両辺を微分したものとして、下から２つ目の式、

```math
1 - F(x) = -\mathbb{E}(T)f(x)
```

が、いきなり「先生」から出されます。

しかし、これは実際に計算したものでは

```math
-1 + F(x) = -\mathbb{E}(T)f(x)
```

となるもので、書籍中のとは符号が違っています。

---

そしてその次の一番下の等式では、

```math
\frac{\frac{d}{dx}(1-F(x))}{1-F(x)} = -\frac{1}{\mathbb{E}(T)}
```

となっています。ここでこの等式をもとに、下から二番目の式の形にもどしても、

```math
-(1 - F(x)) = \mathbb{E}(T)\frac{d}{dx}(1-F(x))
```

となるので、この面からも書籍での下から二番目の式が違っていることは確認できるでしょう。

---

また17ページの

```math
f(x) = \frac{d}{dx}F(x)
```

から、

```math
\frac{d}{dx}(1-F(x)) = -f(x)
```

となることで、この結果を下から二番目では右辺に適用されていますが、そこから導出する一番下の式で微分型に戻しているので、これを適用する必要もないのではと思いました。

---

そもそも2つ目の等式の左辺の微分は、書籍でのように「先生」が直接だして、それを生徒たちもすぐ受け入れられる単純さでもないと思いました。実際、読むことで理解させるレベルにするなら、

```math
\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ページにある、

```math
\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$を想定すると、

```math
F(x) = 1 - e^{-\lambda x}
```

という分布関数が得られ、この分布が「[指数分布](https://ja.wikipedia.org/wiki/%E6%8C%87%E6%95%B0%E5%88%86%E5%B8%83)」とよばれるものである、といった結末になっています。

---

## 指数分布の性質をプログラムで確認する

最後に、プログラミングの話題となるよう、プログラム実行によってこの**時刻にかかわらず待ち時間期待値が同じ**である性質を確認してみます。ここではJavaScriptを使います。

通常の言語環境にビルトインされている乱数機能は一様分布です。
一様分布を用いてある分布関数に従う乱数値を作るには、**その分布の逆関数へ0から1の範囲の一様乱数の値を渡す**ことで実現できます。分布関数はマイナス無限大から引数値までの範囲への累積確率で、0から1の範囲で結果を返すからです。

指数分布は$F(x) = 1-e^{-\lambda x}$なので、この逆関数は$F^{-1}(y) = \frac{-\log(1-y)}{\lambda}$となります。
よって、指数分布に従うランダムな到着時間を返す関数は、以下のようになります:

```js:in[1]
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$ではなく、その逆数である待ち時間の期待値で与えるようにしてあります。

続いて、指数分布の到着時刻リストの内容を一覧するために、数値配列を対数ヒストグラムでコンソールにプロットするコードも用意しておきます。

```js:in[2]
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`での待ち時間の平均を計算して表示するコードは、以下のようになります。

```js:in[3]
{
  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つのコードを上から順に実行した結果は、以下のようになりました:

```txt:out[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...を掛けた値です。

時刻が後ろに行くほど対象となるサンプル数が減ってしまうので理想から外れていきますが、サンプル数が十分あることを想像すれば、どの時刻でも平均値が期待値に近いものとなることがわかるでしょう。