2021-06-10

分布に従う乱数の作り方

一様分布に従う確率変数 $T \sim U[0,1)$ は、多くのプログラミング言語で容易に表現できる。たとえばPythonならば、random モジュールの random.random() を用いると $I = [0, 1)$ 上の一様乱数を得ることができる。

つまり、確率変数を $T = \text{id}:\,\mathbb{R} \rightarrow I, x \mapsto x$ とするとき、 $x \in \mathbb{R}$ に対して

$\displaystyle \mathbb{P}(T \leqq x) = \int _ {-\infty} ^ {x} {\mathbf{1} _ {[0,1)}(t) \text{d}t}$

とかくことができる。ここに、集合 $A$ に対して $\mathbf{1}_A$ は $A$ の定義関数（特性関数）である：

$\mathbf{1}_{A}:\ \mathbb{R} \rightarrow \{0,1\},\ x \mapsto \begin{cases} 1\ \ \ \text{if}\ \ \ x \in A,\\ 0 \ \ \ \text{otherwise.} \end{cases}$

実は、 $\mathbf{1}_{[0,1)}$ は不連続点の集合が高々可算なので、 $\text{Riemann}$ 可積分である。

ここで、ある分布の分布関数を $F$ とする。たとえば、正規分布 $\text{N}(\mu, \sigma ^ 2)$ の確率密度関数は $f(x) = \dfrac{1}{\sqrt{2\pi}\sigma} e ^ {-(x-\mu) ^ 2 / 2\sigma ^ 2}$ で与えられるから、その分布関数 $F(x)$ は、

$\displaystyle F(x) = \dfrac{1}{\sqrt{2\pi}\sigma} \int _ {-\infty} ^ {x} {e ^ {-(t-\mu) ^ 2 / 2\sigma ^ 2} \text{d}t}$

で表される。これは初等関数で表現できないことが知られている。

いま、 $F$ を狭義単調増加であると仮定すると、逆関数 $F ^ {-1}:\ [0,1] \rightarrow \mathbb{R}$ が存在して $F ^ {-1}$ も狭義単調増加であり、 $T(x) = x \in [0,1)$ であるから、 $T(x) \leqq F(x) \Longleftrightarrow F ^ {-1}(T(x)) \leqq x$ が成立つ。なお、 $F$ は分布関数なので単調増加、つまり単調非減少ではある。しかし、狭義単調増加でない場合には、一般に逆関数が存在しないので、この「 $F$ が狭義単調増加」という仮定は必要である（実は、この仮定がなくても、少し工夫すれば「逆関数」に相当する関数を定められるのだが、本稿では省略する）。

さて、このもとで、 $\mathbb{P}(F ^ {-1}(T) \leqq x)$ を計算してみる。なお、この $\mathbb{P}$ の中にある $F ^ {-1}(T) \leqq x$ は $\{F ^ {-1}(T) \leqq x\}$ すなわち $\{x \in I;\ F ^ {-1}(T(x)) \leqq x\}$ を略した記号である（つまり、「関数 $\leqq$ 実数値」の不等式が $\mathbb{P}$ に与えられているのではない）：

$\displaystyle \mathbb{P}(F^{-1}(T) \leqq x) = \mathbb{P}(T \leqq F(x)) = \int _ {-\infty} ^ {F(x)} {\mathbf{1} _ {[0,1)}(t) \text{d}t}$

ここで、 $0 \leqq F \leqq 1$ だから、この積分はかなり簡明になる：

$\displaystyle \int _ {-\infty} ^ {F(x)} {\mathbf{1} _ {[0,1)}(t) \text{d}t} = \int _ {0} ^ {F(x)}\text{d}t = F(x).$

したがって、 $\mathbb{P}(F ^ {-1}(T) \leqq x) = F(x)$ が成立つ。これは、確率変数 $F ^ {-1} \circ T$ が、分布関数が $F$ で表される分布に従うことを示している。なお、ここでは「分布関数と分布が1対1に対応する」という定理を用いており、この事実から、分布関数と分布は同一視されることが多い。つまり、赤字の部分は、数式では

$F ^ {-1} \circ T \sim F$

とかける。すなわち、一様乱数 $x \in [0,1)$ を $F ^ {-1}$ で歪めた $F ^ {-1}(x)$ は、分布 $F$ に従う。

ただし、この方法の難点は、たとえば正規分布であれば、複雑な分布関数の逆関数が必要な点である。

例

Pythonを用いて、正規分布 $\text{N}(\mu, \sigma ^ 2)$ に従う乱数を生成しよう。なお、特殊な記法やライブラリを用いないので、大抵のプログラミング言語で同様のプログラムを書くことができると思う。

先に述べた通り、分布関数は

$\displaystyle F(x) = \dfrac{1}{\sqrt{2\pi}\sigma} \int _ {-\infty} ^ {x} {e ^ {-(t-\mu) ^ 2 / 2\sigma ^ 2} \text{d}t}$

であるが、これの逆関数を得るのは明らかに面倒そうで、特殊なライブラリを用いる必要があったりする。たとえばPythonでは scipy.stats の norm.ppf 関数がこれにあたるが、個人的に手元の環境で pip が正常に動かないので、scipy がインストールできずに詰まった。直そうとは思うのだが、沼に嵌りそうで、比較的忙しくPCを多用するいま、環境を無闇に変えたくなくて、あまりやりたくない。こういうことが gcc についても前にあったのだが、直そうと思って色んなところをいじくり回していたら gcc が #include <stdio.h> を認識できなくなってしまって、それ以降ちょっとトラウマである。

話が逸れた。統計に強いPythonですらこれなのだから、他のプログラミング言語ではそもそも提供されているかも怪しい。いちいち探すのも面倒だから、普遍的な実装を一つ知っておいて、それで代用できるようになっておけばいい。

さて、その問題である逆関数であるが、数値解析、俗に言うNewton法（Newton-Raphson法）を用いる。こう書くとものすごく難しそうだが、その理論はさておいて実装は軽い：

def nr(f, y, df, x0, e):
    x = x0
    while (abs(f(x) - y) >= e):
        x = x - (f(x) - y) / df(x)
    return x

この関数は、簡単に言うと $f(x) = y$ の根 $x ^ *$ の近似解を求める。適当な初期値 $x _ 0$ から始めて、第 $n$ 近似値を $G(x) = f(x) - y$ の $x = x _ {n-1}$ における接線の $x$ 切片で与える：

$x _ n = x _ {n-1} - \dfrac{G(x _ {n-1})}{G ^ {\prime}(x _ {n-1})} = x _ {n-1} - \dfrac{f(x _ {n-1}) - y}{f ^ {\prime}(x _ {n-1})}$

$e$ を $0$ に近い正値とし、この第 $n$ 近似値 $x _ n$ が $|f(x _ n) - y| \lt e$ となった場合に近似を打ち切る、といった手法である。

ところで、このように、一律に閾値を定めて近似精度の尺度とするのはよくあることなのだが、実はあまりよくないことが知られている。 $f$ によっては、 $|f(x _ n) - y| \lt e$ だからといって $x _ n$ が $0$ に近いとは限らないからである。たとえば $G$ が $x$ 軸付近で横這いだった場合を考えるとよい。しかし、今回は簡単のためにこの方法を採る。

また、分布関数

$\displaystyle F(x) = \dfrac{1}{\sqrt{2\pi}\sigma} \int _ {-\infty} ^ {x} {e ^ {-(t-\mu) ^ 2 / 2\sigma ^ 2} \text{d}t}$

は、誤差補関数 $\text{erfc}$ を用いて $F(x) = \dfrac{1}{2}\text{erfc}\left(-\dfrac{(x-\mu)}{\sqrt{2}\sigma}\right)$ とかける。誤差関数や誤差補関数は多くの言語で市民権を得ているように思う。また $F ^ \prime$ は、積分の微分だから簡単で、

$F ^ {\prime}(x) = \dfrac{1}{\sqrt{2\pi}\sigma}\text{exp}\left(-\dfrac{(x-\mu) ^ 2}{2\sigma ^ 2}\right)$

となる。以上が分かれば、プログラムは次のようにかける：

# coding: utf-8

import math
from random import random
import functools

def nr(f, y, df, x0, e):
    x = x0
    while (abs(f(x) - y) >= e):
        x = x - (f(x) - y) / df(x)
    return x
def f(e, v, x):
    return math.erfc(-(x-e) / ((2*v)**(1/2))) / 2
def df(e, v, x):
    return math.exp(-(x-e)**2 / (2*v)) / ((2*math.pi*v)**(1/2))

def main():
    M    = int(input("size:\t"))
    mean = float(input("mean:\t"))
    var  = float(input("var:\t"))
    g  = functools.partial(f, mean, var);
    dg = functools.partial(df, mean, var);

    print("Now sampling...")
    sample = [nr(g, random(), dg, mean, 10**(-5)) for _ in range(M)]
    exp = sum(sample) / M
    var = sum([(x - exp)**2 for x in sample]) / M
    
    print("Exp:\t", exp)
    print("Var:\t", var)

if __name__ == '__main__':
    main()

functools.partial というのを使っているが、これは多変数関数に特定の引数を指定した、新しい関数を作成するへルパ（helper）関数である。たとえば、 $f(x, y)$ に $x = a$ のみを代入すると、 $y$ に関する関数 $f(a, y)$ ができる。この関数 $f(a, y)$ を表現するのが functools.partial(f, a) である。これは便法のためであり、使わなくとも同等のコードは作れる。

実行結果は次の通り：

$ python3 nr.py
size:   10000
mean:   1
var:    25
Now sampling...
Exp:     0.9657192785266343
Var:     24.33909416336915

指定した平均および分散を有するように、ランダムな数を生成できていることが分かる。

2021-03-26

加比の理の条件

加比の理を、

$a, b, c, d \in \mathbb{R}$ について, $a, c, a + c \neq 0$ かつ $\dfrac{b}{a} \lt \dfrac{d}{c}$ ならば $\dfrac{b}{a} \lt \dfrac{b + d}{a + c} \lt \dfrac{d}{c}$ が成立つ

とする人がいますが、例えば $\dfrac{2}{3} \lt 1 = \dfrac{-1}{-1}$ だから $\dfrac{2}{3} \lt \dfrac{2-1}{3-1} = \dfrac{1}{2} \lt 1$ となるので誤りです。

いま、

$\dfrac{b}{a} \lt \dfrac{d}{c}$ $\Longleftrightarrow$ $\dfrac{d}{c} - \dfrac{b}{a} = \dfrac{ad - bc}{ac} \gt 0$ （＊）

が成立するので、

$\dfrac{b}{a} \lt \dfrac{b + d}{a + c}$ $\Longleftrightarrow$ $\dfrac{b + d}{a + c} - \dfrac{b}{a} = \dfrac{ad - bc}{a(a + c)} \gt 0$ （＊＊）

であることから、

① $a \gt 0,\ c \gt 0$ なら（このときは $a + c \gt 0$ ）、（＊）から $\dfrac{ad - bc}{a} \gt 0$ であるので（＊＊）の右辺が成立して○。

② $a \lt 0,\ c \lt 0$ なら（このときは $a + c \lt 0$ ）、同様に考えて○。

③ $a \gt 0,\ c \lt 0$ なら、 $a + c \gt 0$ のときは○で、 $a + c \lt 0$ のときは×。

④ $a \lt 0,\ c \gt 0$ なら、 $a + c \gt 0$ のときは×で、 $a + c \lt 0$ のときは○。

また、同様にして、

$\dfrac{d}{c} \gt \dfrac{b + d}{a + c}$ $\Longleftrightarrow$ $\dfrac{d}{c} - \dfrac{b + d}{a + c} = \dfrac{ad - bc}{c(a + c)} \gt 0$ （＊＊＊）

であることから、

① $a \gt 0,\ c \gt 0$ なら（このときは $a + c \gt 0$ ）、○。

② $a \lt 0,\ c \lt 0$ なら（このときは $a + c \lt 0$ ）、○。

③ $a \gt 0,\ c \lt 0$ なら、 $a + c \gt 0$ のときは×で、 $a + c \lt 0$ のときは○。

④ $a \lt 0,\ c \gt 0$ なら、 $a + c \gt 0$ のときは○で、 $a + c \lt 0$ のときは×。

以上から、（＊＊）かつ（＊＊＊）が成立つための十分条件として、「 $a \gt 0,\ c \gt 0$ または $a \lt 0,\ c \lt 0$ 」 $\Longleftrightarrow$ $ac > 0$ であればよいので、加比の理は、

$a, b, c, d \in \mathbb{R}$ について, $a, c, a + c \neq 0$ かつ $ac > 0$ かつ $\dfrac{b}{a} \lt \dfrac{d}{c}$ ならば $\dfrac{b}{a} \lt \dfrac{b + d}{a + c} \lt \dfrac{d}{c}$ が成立つ

とかけます。

2021-02-19

包除の定理

〔参考文献：Williram Feller（河田龍夫ら訳）『確率論とその応用Ⅰ（上巻）』紀伊国屋書店, 2001.〕

（離散）標本集合を $\Omega$ とし、その部分集合である $n$ 個の事象 $A _ 1,\ A _ 2,\ \cdots,\ A _ n$ を考える。「ちょうど $m$ 個の事象のみが起こる」という事象 $A$ の起こる確率は、

$\displaystyle S _ k = \sum _ {\sigma \in \{ \sigma \in \mathfrak{S} _ k\ |\ \sigma(1) \lt \cdots \lt \sigma(k) \} } {\mbox{Pr}(A _ {\sigma(1)} \cap \cdots \cap A _ {\sigma(k)})}$

とおくとき（ $\mathfrak{S} _ n$ は $n$ 次の置換全体の集合）、

$P _ {[m]} = S _ m - \binom{m+1}{m} S _ {m+1} + \binom{m+2}{m} S _ {m+2} - \cdots \pm \binom{n}{m} S _ {n}$

で求めることができる。長い式で少し嫌になる（特にΣで和をとる範囲が面倒そうに思える）のだが、実用の幅が広い定理で（冒頭に示した参考文献で、色々の応用を扱っている）、とても重要だと思ったので、メモする。

証明. 次の式で、指示関数（indicator） という関数：

$\boldsymbol{1} _ {A} (x) = 1 \ (x \in A),\ 0 \ (x \notin A)$

を定める（ $\boldsymbol{1} _ A$ は $\chi _ A$ などとも書かれる）。つまり $\boldsymbol{1} _ A$ は、 $x$ を入力とし、 $x \in A$ ならば $1$ を, $x \notin A$ ならば $0$ を返すような関数である。この指示関数を用いると、

$\displaystyle P _ {[m]} = \sum _ {x \in A} {\mbox{Pr}(\{x\})} = \sum _ {x \in \Omega} {\mbox{Pr}(\{x\})} \boldsymbol{1} _ {A}(x),$

$\displaystyle S _ k = \sum _ {\sigma \in \{ \sigma \in \mathfrak{S} _ k\ |\ \sigma(1) \lt \cdots \lt \sigma(k) \} } { \sum _ {x \in \Omega} \mbox{Pr}(\{x\}) \boldsymbol{1} _ {A _ {\sigma(1)} \cap \cdots \cap A _ {\sigma(k)}}(x)}$ $\displaystyle = \sum _ {x \in \Omega} \mbox{Pr}(\{x\}) \sum _ {\sigma \in \{ \sigma \in \mathfrak{S} _ k\ |\ \sigma(1) \lt \cdots \lt \sigma(k) \} } \boldsymbol{1} _ {A _ {\sigma(1)} \cap \cdots \cap A _ {\sigma(k)}}(x)$

と表すことができる。記号の便宜で、

$\displaystyle s _ k(x) = \sum _ {\sigma \in \{ \sigma \in \mathfrak{S} _ k\ |\ \sigma(1) \lt \cdots \lt \sigma(k) \} } \boldsymbol{1} _ {A _ {\sigma(1)} \cap \cdots \cap A _ {\sigma(k)}}(x)$

とおくことにしよう。よって $S _ k = \sum _ {x \in \Omega} \mbox{Pr}(\{x\}) s _ k(x)$ となる。以上のことから、示したい式は

$\sum _ {x \in \Omega} {\mbox{Pr}(\{x\})} \boldsymbol{1} _ {A}(x)$

$= \sum _ {x \in \Omega} \mbox{Pr}(\{x\}) (s _ m(x) - \binom{m+1}{m} s _ {m+1}(x)$ $\hspace{1in} + \binom{m+2}{m} s _ {m+2}(x) - \cdots \pm \binom{n}{m} s _ {n}(x))$

と書くことができるので、任意の $x \in \Omega$ について、

$\boldsymbol{1} _ A(x)$ $= s _ m(x) - \binom{m+1}{m} s _ {m+1}(x) + \binom{m+2}{m} s _ {m+2}(x) - \cdots \pm \binom{n}{m} s _ {n}(x)$

を示すことができれば証明が終わる（もはやこれは確率論というよりも集合論であることに注意されたい。この定理の本質はこの式に他ならない。 なお、例えば $\Omega$ を有限集合と仮定して、 $T \subseteqq \Omega$ について $\mbox{n}(T)$ で $T$ に含まれる元の個数を表すことにすれば、 $\mbox{Pr}()$ を $\mbox{n}()$ で置き換えるだけで上の議論をそのまま使って $\mbox{n}()$ に関する同様の定理が証明できる）。

$x \in A$ の場合。このとき、 $x$ はちょうど $m$ 個の事象 $A _ {i _ 1},\ \cdots,\ A _ {i _ m}$ のみに含まれるのだったから、 $s _ {m+1}(x) = \cdots = s _ {n}(x) = 0$ となるので、 $s _ {m}(x) = 1$ を示せばよい。 $i _ 1,\ \cdots,\ i _ m$ を昇順（小さい順）に並び替えて $j _ 1 \lt \cdots \lt j _ m$ が一意に定まるので、 $s _ {m}(x)$ のΣのうちで、 $\sigma(1) = j _ 1 ,\ \cdots,\ \sigma(m) = j _ m$ なるただひとつの $\sigma$ についての項のみが残り、この項の値は $1$ である。ゆえに $s _ m(x) = 1$ となる。

$x \notin A$ の場合。 $x$ は $A _ 1,\ \cdots,\ A _ n$ のうちで、ちょうど $r$ 個 $(r \neq m)$ の事象に含まれるとする。 $r \lt m$ ならば $s _ m(x) = \cdots = s _ {n}(x) = 0$ となるから成立するので、以下 $r \gt m$ とする。なお、やはり $s _ {r + 1}(x) = \cdots = s _ {n}(x) = 0$ は成立つので、示すべくは、

$0 = s _ {m}(x) - \binom{m+1}{m} s _ {m+1}(x) + \cdots \pm \binom{r}{m} s _ {r}(x)$

である。

$x$ は $r\ ( \gt m)$ 個の事象に属するから、 $m \leqq c \leqq r$ について、 $s _ {c}(x)$ のΣのうちで、残る項は $\binom{r}{c}$ 個だけあり（これがポイントである。 $\sigma(1) \lt \cdots \lt \sigma(c)$ という条件が附されているので、 $r$ 個から $c$ 個を選ぶ組合せの総数と同じだけ項が残る）、その値はすべて $1$ であるから、 $s _ {c}(x) = \binom{r}{c}$ が成立つ。したがって、示す式は、

$0 = \binom{r}{m} - \binom{m+1}{m}\binom{r}{m+1} + \cdots \pm \binom{r}{m}\binom{r}{r}$

とかける。ここで $\binom{m+y}{m} \binom{r}{m+y} = \binom{r}{m} \binom{r-m}{y}$ が成立するので、

$0 = \binom{r}{m} \left\{\binom{r-m}{0} - \binom{r-m}{1} + \cdots \pm \binom{r-m}{r-m}\right\}$

と変形できるが、この右辺の $\{\}$ の中は $(1+(-1)) ^ {r-m}$ の二項展開だから（右辺） $=0$ が得られる。これで証明が終わった。終

この定理の効力を見るために、参考文献から一例を借用しよう。

$1, 2, \cdots, n$ と番号の書かれたカードが二組あって、一方を番号が昇順になるように、横一列に机に並べる。そのすぐ下に、他方をシャッフルして横一列に並べる。結果として机には二列の $n$ 枚のカードが並び、上は整列されていて、下は勝手に並んでいる、という状況になる。 $n= 5$ ならば、例えば

$1\ \ 2\ \ 3\ \ 4\ \ 5$

$3\ \ 2\ \ 1\ \ 5\ \ 4$

という風に並んでいるのである。このとき、この二列を縦にみて、上下で数字が一致しているか否かを考えよう。この例では、 $2$ のみが一致している。

さて、下列の並び方は全部で $n !$ 通りある。そこで、どの並び方も同じ確率 $1 / n !$ で起こると仮定するとき、ちょうど $m$ 個の数字が一致する確率はどれくらいだろうか？

下列の並び方を、数字をカンマ区切りで横に並べて括弧でまとめることで表現しよう。先の例なら $(3,2,1,5,4)$ と表すことができる。逆に、 $n$ 個の数字によるこの表記から、ある一つの下列の並べ方が対応するので、下列の並び方のすべてをこのような表記で統一的に表すことができる。従って、直積を用いて標本空間 $\Omega$ は $\Omega = \{1, 2, \cdots, n\} ^ {n}$ とかける。

事象 $A _ k$ を「番号 $k$ が一致する」ことによって定める $(1 \leqq k \leqq n)$ 。厳密に言えば $A _ k = \{(i _ 1, i _ 2, \cdots, i _ {k-1}, k, i _ {k+1}, \cdots, i _ {n})\ |\ 1 \leqq i _ 1, i _ 2, \cdots, i _ n \leqq n\}$ とかくことができる。含まれる元の個数は $(n-1) !$ である。求めたい確率は、事象 $A$ を「 $A _ 1, \cdots, A _ n$ のうちでちょうど $m$ 個のみの事象が起こる」ことと定義して $P _ {[m]} = \mbox{Pr}(A)$ とかけるので、定理の記号で、

$\displaystyle S _ k = \binom{n}{k} \dfrac{(n-k) !}{n !} = \dfrac{1}{k !}$

であるから、

$P _ {[m]} = S _ m - \binom{m+1}{m} S _ {m+1} + \binom{m+2}{m} S _ {m+2} - \cdots \pm \binom{n}{m} S _ {n}$

$= \dfrac{1}{m !} - \dfrac{(m+1) !}{m! 1 !} \dfrac{1}{(m+1) !} + \dfrac{(m+2) !}{m! 2 !} \dfrac{1}{(m+2)!} - \cdots \pm \dfrac{n !}{m!(n-m)!}\dfrac{1}{n!}$

$= \dfrac{1}{m !} \left(\dfrac{1}{0!} - \dfrac{1}{1!} + \dfrac{1}{2!} - \cdots \pm \dfrac{1}{n!}\right)$

$\displaystyle = \dfrac{1}{m!} \sum _ {k = 0} ^ {n} {\dfrac{(-1) ^ k}{k !}}$

と求めることができる。特に $n$ が十分大ならば、 $\displaystyle \sum _ {k = 0} ^ {\infty} {\dfrac{(-1) ^ k}{k !}} = e ^ {-1}$ であることから、

$P _ {[m]} \sim \dfrac{e ^ {-1}}{m !}$

と近似できる。さらに $m = 0$ とすれば、

$P _ {[0]} \sim e ^ {-1} = 0.36787...$

であるが、これは「十分多い $n$ 枚のカードのもとで、ひとつも数字が一致しない並び方となる確率」を表している。

$\sigma \in \mathfrak{S} _ n$ は、すべての $1 \leqq i \leqq n$ について $\sigma(i) \neq i$ となるとき撹拌置換と呼ばれる。この言葉を使えば、 $n$ 次の撹拌置換を得る確率は概ね $37 \%$ 程度である。

2021-02-18

三囚人問題

有名問題。

ある刑務所長は3人の囚人の中からランダムに1人を選んで解放し, 残りの2人を処刑する. 看守は誰が解放されるか知っているが, どの囚人にも彼が開放されるかどうか教えることは禁止されている. 囚人を $X, Y, Z$ と呼ぶことにしよう. $X$ は看守に $Y$ か $Z$ のどちらが処刑されるか尋ね, 私は $Y$ か $Z$ のどちらかが処刑されることは知っているので, 看守がどちらが処刑されるかを教えても, 私の状態について情報を漏らしたことにはならないと主張した. 看守は $X$ に $Y$ が処刑されることを伝えた. $X$ は, 彼か $Z$ のどちらかが解放されるので, 彼が解放される確率は今や $1/2$ であり, 以前より幸せに感じている. 彼は正しいか, それとも彼のチャンスはまだ $1/3$ のままか？説明せよ.

〔出典：T. H. Cormen, C. E. Leiserson, R. L. Rivest, C. Stein: Introduction to Algorithm (3rd Edition), MIT Press, 2009. （浅野哲夫, 岩野和生, 梅尾博司, 山下雅史, 和田幸一（訳）『アルゴリズムイントロダクション第3版』近代科学社, 2013）〕

問題について

モンティホール問題と同様に、一見してどっちが正しいのか分からなくなるような問である。

＊

最初に正しい考え方を示して、その後に誤った考え方を示そう。

確率を議論するためには、起こり得るすべての結果の全体の集合である、標本空間 $\Omega$ を定めねばならない。起こり得るすべての結果は、それぞれ（解放される囚人, 看守の情報）という二つ組で表現できるので、

$\Omega = \{(X,Y), (X,Z), (Y,Z), (Z,Y)\}$

と表せる。いま、解放される囚人は同様に確からしいから、

$\mbox{Pr}(\{(X,Y), (X,Z)\}) = \mbox{Pr}(\{(Y,Z)\}) = \mbox{Pr}(\{(Z,Y)\})$

が成立つ。 $\mbox{Pr}(\Omega) = 1$ であるから、これら三つの確率は $1/3$ で等しい。さらに、看守が無作為の選択をすることに注意する。解放される囚人が $X$ であった場合は $Y$ も $Z$ も処刑されるので、看守は $X$ に伝える囚人として、 $Y$ と $Z$ のうちどちらか一方を選択しなければならないが、問題の記述からこれは無作為であると仮定せざるを得ない。したがって、

$\mbox{Pr}(\{(X,Y)\}) = \mbox{Pr}(\{(X,Z)\})$

が成立つ。ここで、 $\mbox{Pr}(\{(X,Y), (X,Z)\}) = 1/3$ であったから、これら二つの確率は $1/6$ で等しい。以上をまとめると、次のようになる。

$\mbox{Pr}(\{(X,Y)\}) = \mbox{Pr}(\{(X,Z)\}) = 1/6,$ $\mbox{Pr}(\{(Y,Z)\}) = \mbox{Pr}(\{(Z,Y)\}) = 1/3.$

求める確率は、「 $A:$ 看守が $Y$ を伝える」もとで、「 $B:$ $X$ が解放される」ような条件付確率 $\mbox{Pr}(B\ |\ A) = \mbox{Pr}(A \cap B) / \mbox{Pr}(A)$ である。 $A = \{(X,Y), (Z,Y)\},$ $B = \{(X,Y), (X,Z)\}$ であるので、

$\mbox{Pr}(B\ |\ A) = \mbox{Pr}(\{(X,Y)\}) / \mbox{Pr}(\{(X,Y), (Z,Y)\}) = (1/6)/(1/6+1/3) = \boldsymbol{1/3}.$

つまり $X$ の解放される確率は看守の助言によって変化しない。さて、次のように考えるのはどうだろうか？

我々は「起こり得るすべての結果は、それぞれ（解放される囚人, 看守の情報）という二つ組で表現できる」として話を進めたが、結果として看守は「 $Y$ が処刑される」と言うのだから、結果を表すのに看守の情報はもはや必要なく、標本空間 $\Omega$ は $\Omega = \{X, Z\}$ と考えてよい。解放される囚人は無作為に決定されるのだから、 $\mbox{Pr}(\{X\}) = \mbox{Pr}(\{Z\})$ であり、 $\mbox{Pr}(\Omega) = 1$ であるので $\mbox{Pr}(\{X\}) = 1/2$ となるから、 $X$ の考えは正しい。

この考え方は誤っている。看守が「 $Y$ が処刑される」と助言できるのは、解放される囚人を無作為に選んだ後にのみ可能であることに注意すると、この誤答では因果関係が逆転していることに気づくだろう。看守の助言は、解放される囚人の無作為な選択に依存して定まる。確率を定義するためには、それ以前に標本空間が与えられていることを前提とするので、これでは最初の「無作為に解放囚人を選ぶ」とき以前に標本空間が与えられていないことになって不合理を生ずる。つまり、「標本空間を定める」→「無作為に解放囚人を選ぶ」→「看守から助言を受ける」が正しい依存関係であるのに、「無作為に解放囚人を選ぶ」→「看守から助言を受ける」→「標本空間を定める」としてしまっているので意味がない。標本空間は過程の初めに定まっていなければならない。 これが確率の定義であり、このことを忘れると結果を誤りがちである。

2021-02-17

モンティホール問題

有名問題。

あなたはあるクイズショウの参加者である. 3つあるカーテンのうちの1つの後ろに賞品が隠されていて, 正しいカーテンを選べばこの賞品を獲得できる. あなたが1つのカーテンを選んだ後, そのカーテンを引き上げる前に司会者は残りのカーテンの中から1つを引き上げ, 空であることを示して, あなたに決心を変えるかどうか尋ねる. あなたが決心を変えたとき, あなたが賞品を得るチャンスは変わるか？

問題について

多くの人が、「最終的に二択になるんだから単純に半々の確率では？」と思うようで、実際自分もそう思ったが、違うらしく、正解は「選び直せば2/3、選び直さなければ1/3で当てることができる」という、直観にいくらか反したものである。具体的に計算すれば確かにその通りになる。詳しくは↓

mathtrain.jp

「ああ確かに」と納得してしまえばそれまでだが、もう少し考えてみた。正しい答は知ったが、なぜ間違えたのかという本質的な問題がしばらく解消できなかった。つまり、やはり直観による答も正しいように思えて、直観の何が正しく何が間違っているのかが判然としなかった。今日ずっと考えて、自分なりに整理できたような気がするので、まとめておきたい。

＊

確率を定めるためには、起こり得るすべての結果を含む集合である標本空間 $\Omega$ を表現しなければならない（これは当たり前のように聞こえるが、とても重要なことである）。このゲームの進み方は

①最初に無作為に選んだカーテンの番号

②司会者から教わった外れのカーテンの番号

③最後に選んだカーテンの番号

という三つの数の組合せで表現することができる。カーテンの番号を順に $0, 1, 2$ 番とする。例えば $(0,1,0)$ は、①最初に0番のカーテンを選び、②司会者から1番が外れだと教わり、③決心を変えずに0番を選んだ－という進行を表現するものとしよう。いま、このゲームで言う「正しいカーテン」が0番であると固定すると、進行の全体の集合 $\Omega$ は次のようにかける：

$\Omega = \{(0,1,0), (0,1,2), (0,2,0), (0,2,1), (1,2,0), (1,2,1), (2,1,0), (2,1,2)\}.$

ここが誤解だが、これらすべての結果が同確率で起こる（つまり、「同様に確からしい」）わけではないことに注意したい。「同様に確からしい」という概念は、無作為に選んだ結果に仮定できる条件である。我々は①で、最初の番号を無作為に選んでいるので、

$\mbox{Pr}(\{(0,1,0), (0,1,2), (0,2,0), (0,2,1)\})$ $= \mbox{Pr}(\{(1,2,0), (1,2,1)\}) = \mbox{Pr}(\{(2,1,0), (2,1,2)\})$

を結論づけることはできる。しかし、

$\mbox{Pr}(\{(0,1,0)\}) = \mbox{Pr}(\{(0,1,2)\}) = \mbox{Pr}(\{(0,2,0)\})$ $= \mbox{Pr}(\{(0,2,1)\}) = \mbox{Pr}(\{(1,2,0)\}) = \mbox{Pr}(\{(1,2,1)\})$ $= \mbox{Pr}(\{(2,1,0)\}) = \mbox{Pr}(\{(2,1,2)\})$

というように、すべての結果が同確率で起こると仮定するには根拠がない。とはいえ、「根拠がない」と一蹴して終わるのも味気ないので、実際にこの仮定の下で話を進めて、間違っていることを見てみよう。確率の性質によって、 $\mbox{Pr}(\Omega) = 1$ であり、部分集合 $A, B \subseteqq \Omega$ について $A \cap B = \emptyset$ ならば $\mbox{Pr}(A\cup B) = \mbox{Pr}(A) + \mbox{Pr}(B)$ が成立つことを用いると、

$\mbox{Pr}(\Omega) =$ $\mbox{Pr}(\{(0,1,0)\}) + \mbox{Pr}(\{(0,1,2)\}) + \mbox{Pr}(\{(0,2,0)\})$ $+ \mbox{Pr}(\{(0,2,1)\}) + \mbox{Pr}(\{(1,2,0)\}) + \mbox{Pr}(\{(1,2,1)\})$ $+ \mbox{Pr}(\{(2,1,0)\}) + \mbox{Pr}(\{(2,1,2)\})$ $= 1$

であって、すべての結果が同確率で起こるという仮定から、すべての結果は確率 $1/8$ で起こることが分かる。すると、番号を変えて当たる確率は

$\mbox{Pr}(\{(1,2,0), (2,1,0)\}) = 1/8 + 1/8 = 1/4$

であり、番号を変えずに当たる確率も

$\mbox{Pr}(\{(0,1,0), (0,2,0)\}) = 1/8 + 1/8 = 1/4$

となって、番号を変えようが変えまいが同じ確率になって、正しい答が得られていないのが分かる。

本題に戻る。先のことから、

$\mbox{Pr}(\{(0,1,0), (0,1,2), (0,2,0), (0,2,1)\})$ $= \mbox{Pr}(\{(1,2,0), (1,2,1)\}) = \mbox{Pr}(\{(2,1,0), (2,1,2)\})$

は成立つのだった。 $\mbox{Pr}(\Omega) = 1$ だから、これら三つの確率は $1/3$ で等しいことが分かる。ここで、司会者も無作為の選択をしていることに注意する。挑戦者が見事①で正解 $0$ を選んだ場合、司会者が空ける外れのカーテンの番号の候補は二つあり、問題の仮定から、司会者はこの二つから空けるカーテンの番号を同確率で選ぶと考えざるを得ない。したがって、①で0が選ばれた場合に、②で選ばれる進行は同様に確からしい。ゆえに、

$\mbox{Pr}(\{(0,1,0), (0,1,2)\}) = \mbox{Pr}(\{(0,2,0), (0,2,1)\})$

が成立つので、 $\mbox{Pr}(\{(0,1,0), (0,1,2), (0,2,0), (0,2,1)\}) = 1/3$ であることから、これら二つの確率は $1/6$ で等しい。得られた結果を整理すると、次のようになる。

$\mbox{Pr}(\{(0,1,0), (0,1,2)\}) = 1/6,$

$\mbox{Pr}(\{(0,2,0), (0,2,1)\}) = 1/6,$

$\mbox{Pr}(\{(1,2,0), (1,2,1)\}) = 1/3,$

$\mbox{Pr}(\{(2,1,0), (2,1,2)\}) = 1/3.$

番号を変える場合と変えない場合で、当たる確率と外れる確率をそれぞれ求めよう。

番号を変える場合。 $\mbox{Pr}(\{(0,1,0)\}) = \mbox{Pr}(\{(0,2,0)\})$ $= \mbox{Pr}(\{(1,2,1)\}) = \mbox{Pr}(\{(2,1,2)\}) = 0$ が得られるから、

$\mbox{Pr}(\{(0,1,2)\}) = 1/6,$ $\mbox{Pr}(\{(0,2,1)\}) = 1/6,$

$\mbox{Pr}(\{(1,2,0)\}) = 1/3,$ $\mbox{Pr}(\{(2,1,0)\}) = 1/3.$

従って、「変えて当たる確率」は $\mbox{Pr}(\{(1,2,0), (2,1,0)\}) = \boldsymbol{2/3}$ で、「変えて外れる確率」は $\mbox{Pr}(\{(0,1,2), (0,2,1)\}) = 1/3$ となる。

番号を変えない場合。 $\mbox{Pr}(\{(0,1,2)\}) = \mbox{Pr}(\{(0,2,1)\})$ $= \mbox{Pr}(\{(1,2,0)\}) = \mbox{Pr}(\{(2,1,0)\}) = 0$ が得られるから、

$\mbox{Pr}(\{(0,1,0)\}) = 1/6,$ $\mbox{Pr}(\{(0,2,0)\}) = 1/6,$

$\mbox{Pr}(\{(1,2,1)\}) = 1/3,$ $\mbox{Pr}(\{(2,1,2)\}) = 1/3.$

従って、「変えないで当たる確率」は $\mbox{Pr}(\{(0,1,0), (0,2,0)\}) = \boldsymbol{1/3}$ で、「変えないで外れる確率」は $\mbox{Pr}(\{(1,2,1), (2,1,2)\}) = 2/3$ となる。

司会者の闖入により、最後の③で「決心を変える」「決心を変えない」という二択を迫られるので、「 $1/2$ では？」と思ってしまいがちだが、これは考えている標本空間を取り違えてしまっていることに起因する誤りである。詳しくは三囚人問題を参照のこと。

2021-02-12

ceil(x/(ab))=ceil(ceil(x/a)/b)

面白い証明を見つけたのでメモします。出典は↓

math.stackexchange.com

任意の $x \in \mathbb{R}$ と整数 $a \gt 0$ , $b \gt 0$ について,

$\left\lceil \dfrac{\lceil x / a \rceil}{b} \right\rceil = \left\lceil \dfrac{x}{ab}\right\rceil$

を示せ.

証明. まず, 任意の整数 $n$ と実数 $x$ について, 同値 $n \geqq x \Longleftrightarrow n \geqq \lceil x \rceil$ を示す. $x \leqq \lceil x \rceil$ だから $\Longleftarrow$ は明らか. 逆に $n \geqq x$ ならば, ceil関数の定義から $S =\{m \in \mathbb{Z} ;\ x \leqq m \}$ として $\lceil x \rceil = \min S$ かつ $n \in S$ であるから $n \geqq \lceil x \rceil$ .

したがって, 任意の整数 $k \in \mathbb{Z}$ について, 次の同値が成立つ：

$k \geqq \lceil \lceil x/a \rceil /b \rceil$ $\Longleftrightarrow$ $k \geqq \lceil x /a \rceil / b$ $\Longleftrightarrow$ $bk \geqq \lceil x / a \rceil$ $\Longleftrightarrow$ $bk \geqq x/a$ $\Longleftrightarrow$ $k \geqq x / (ab)$ $\Longleftrightarrow$ $k \geqq \lceil x / (ab) \rceil$ .

特に $\lceil \lceil x/a \rceil /b \rceil \geqq \lceil x / (ab) \rceil$ と $\lceil x / (ab) \rceil \geqq \lceil \lceil x/a \rceil /b \rceil$ が成立つので, 結論を得る. 終

2021-02-01

極限の問題

（問題） $\displaystyle \lim _ {n\to\infty} \left(1-\dfrac{1}{n}\right) ^ {n ^ 2} e ^ {n}$ を求めよ.

（コメント）ついさっき締切られたレポートの問題を解くのに必要になった問です。僕なんかはこれ「 $\left\{\left(1 - \dfrac{1}{n}\right) ^ {-n}\right\} ^ {-n} e ^ {n}$ で、 $n \to \infty$ のとき $e ^ {-n} \times e ^ {n}$ みたいなものだから $1$ だ」と高を括ってしまいましたが、そうじゃないんですね。答を見れば分かりますが、 $\displaystyle \left(1-\dfrac{1}{n}\right) ^ {n ^ 2}$ の方が $e ^ {-n}$ よりも減りが大きいぽいです。同じような問題をやさ理かハイ理で見かけて、同じように間違えた気がします。今ならちゃんと説明できるので、戒めとして記します。

（解答）あとで $t = 1/n$ とおけばよいので, まず $(1-t) ^ {1/{t ^ 2}}$ の $t \rightarrow 0$ における極限を求める. $\log$ をとれば $\dfrac{\log{(1-t)}}{t ^ 2}$ となるが, Taylorの公式から,

（ $\displaystyle \log{(1+x)} = \sum _ {n = 1} ^ {\infty} {\dfrac{(-1) ^ {n-1}}{n}x ^ n}$ であったから）

$\log{(1-t)} = -t -\dfrac{t ^ 2}{2} +o(t ^ 2) \hspace{0.1in} (t \to 0)$

であり, $\dfrac{\log{(1-t)}}{t ^ 2} = -\dfrac{1}{t} - \dfrac{1}{2} + o(1). \hspace{0.1in} (t \to 0)$ よって,

$(1-t) ^ {1/{t ^ 2}} = e ^ {\frac{\log{(1-t)}}{t ^ 2}} = e ^ {-1/t-1/2+o(1)} \Longleftrightarrow (1-t) ^ {1/{t ^ 2}} e ^ {1/t} = e ^ {-1/2 + o(1)}. \hspace{0.1in} (t \to 0)$