而して

ノートとかメモとか。

分布に従う乱数の作り方

一様分布に従う確率変数  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.statsnorm.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

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