プログラマの為の
数学勉強会
第14回

謝辞

この会の企画・会場設備の提供をして頂きました
㈱ ワークスアプリケーションズ様
にこの場をお借りして御礼申し上げます。

この資料について

• http://nineties.github.com/math-seminar に置いてあります。
• SVGに対応したブラウザで見て下さい。主要なブラウザで古いバージョンでなければ大丈夫だと思います。
• 内容の誤り、プログラムのバグは@9_tiesかkoichi.nakamur AT gmail.comまでご連絡下さい。
• サンプルプログラムはPython及びMaximaで記述しています。

擬似乱数の生成法

与えられた確率分布に従う擬似乱数を生成する方法を説明します. 有名な確率分布に従う乱数生成器は標準で用意されていると思いますが, 自前での実装が必要な場面もよくあります.

後の回に説明するモンテカルロ法の基礎となりますので大切です.

前提条件

区間 \([0,1)\) の一様分布(標準一様分布)に従う乱数(標準一様乱数)を生成する, 良質な生成器が存在する事を仮定します. \[ \rho(x) = \left\{\begin{array}{cl} 1 & (0 \leq x < 1) \\ 0 & (\text{上記以外}) \end{array}\right. \]

メルセンヌ・ツイスター法であれば間違いありません. 言語処理系の標準乱数生成器のアルゴリズムは必ず調べるようにしてください.

標準一様分布に従う変数 \(U\) で直接表すことの出来る変数は簡単です.

区間 \([a,b)\) の一様分布に従う変数 \(X\) は \[ X = a + (b-a)U \] と表す事が出来ます. よって,以下の関数で生成する事が出来ます.

from random import random
def uniform(a, b):
    return a + random()*(b-a)

成功確率 \(p\) のベルヌーイ分布に従う変数 \(X\) は \[ X = \left\{\begin{array}{cl} 1 & (U \leq p) \\ 0 & (\text{上記以外}) \end{array}\right. \] と表す事が出来ます.

from random import random
def bernoulli(p):
    if random() <= p:
        return 1
    else:
        return 0

\( X_1,X_2,\ldots,X_n \) が成功確率 \(p\) のベルヌーイ分布に従うならば, 確率変数 \[ X = X_1 + X_2 + \cdots X_n \] は二項分布 \(B(n,p)\) に従うのでした.

よって正規分布に従う乱数を既与とすれば, 以下の関数で \(B(n,p)\)に従う乱数を生成出来ます.

from random import random
def binomial(n,p):
    x = 0
    for i in range(n):
        x += bernoulli(p)
    return x

この方法は \(n\) が大きい場合の計算量が問題となります. その場合は, 中心極限定理によって正規分布で近似する事が出来ます.

また,正規表現に従う変数は以下の定理で変形されます.

【証明】
\(y=ax+b\) と変数変換を行うと, \[ \begin{aligned} \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\} \mathrm{d} x &= \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{-\frac{(y-a\mu-b)^2}{2a^2\sigma^2}\right\} \frac{1}{a}\mathrm{d} y \\ &=\frac{1}{\sqrt{2\pi a^2\sigma^2}}\exp\left\{-\frac{(y-a\mu-b)^2}{2a^2\sigma^2}\right\}\mathrm{d} y \end{aligned} \] □

ベルヌーイ分布は \[ \text{平均}:p,\quad\text{分散}: p(1-p) \] なので, \(X_1,X_2,\ldots,X_n\) がこの分布に従う場合, \[ \overline{X} = \frac{X_1+X_2+\cdots + X_n}{n} \] は \(n\) が大きいとき正規分布 \( N\left(p, \frac{p(1-p)}{n}\right) \) に従います.

つまり, \(X = X_1 + X_2 + \cdots + X_n\) は正規分布 \( N\left(np, np(1-p)\right) \) に従いますので, 以下の関数で二項分布に従う乱数を生成する事が出来ます.

from random import gauss
from math import sqrt
def binomial(n, p):
    return round(gauss(n*p, sqrt(n*p*(1-p))))

逆変換法

続いて, 逆変換法という方法を紹介します.

パラーメータ \(\lambda\) の指数分布の確率密度関数は \[ \rho(x) = \left\{\begin{array}{cl} \lambda e^{-\lambda x} & (x \geq 0) \\ 0 & (x < 0) \end{array}\right. \] なので, この累積分布関数は \(x \geq 0\) のとき \[ \begin{aligned} F(x)&= \int_{-\infty}^x\rho(x)\mathrm{d}x = \int_0^x \lambda e^{-\lambda x}\mathrm{d}x \\ &= [-e^{-\lambda x}]_0^x = \color{yellow}{ 1 - e^{-\lambda x} } \end{aligned} \] となります. \(x < 0\) のときは0です.

以下が \(\lambda = 5\) の指数分布の密度関数のグラフ,

以下がその累積分布関数のグラフです.

\[ F(x) = \left\{\begin{array}{cl} 0 & (x < 1) \\ 1/6 & (1 \leq x < 2) \\ 2/6 & (2 \leq x < 3) \\ 3/6 & (3 \leq x < 5) \\ 4/6 & (4 \leq x < 5) \\ 5/6 & (5 \leq x < 6) \\ 6/6 & (6 \leq x) \\ \end{array}\right. \]

連続的な確率分布に対しては以下が成立します.

【証明】
前半の命題は明らかなので省略.
後半に関して. \(F\) は単調増加関数であるから \[ X \leq x\Leftrightarrow F(X)\leq F(x) \] である. すなわち \[ P(F(X)\leq F(x)) = P(X\leq x) = F(x) \] が成立するので, \(F(x) = x'\quad (0\leq x'\leq 1)\) とおき直せば \[ P(F(X) \leq x') = x'\quad (0 \leq x'\leq 1) \] となる. これは \(F(X)\) が標準一様分布に従う事を示す. □

逆変換法(連続的分布)

既に計算したように, パラーメータ\(\lambda\)の指数分布の累積分布関数は \[ F(x) = 1-e^{-\lambda x}\quad (x \geq 0) \] なので, \[ y = 1-e^{-\lambda x} \Leftrightarrow e^{-\lambda x} = 1-y \Leftrightarrow x = -\frac{1}{\lambda}\ln(1-y) \] より \[ \color{yellow}{ F^{-1}(x) = -\frac{1}{\lambda}\ln (1-x) } \] によって標準一様乱数を変換すればよいです.

\(X\) が標準一様分布に従うならば, \(1-X\) も標準一様分布に従うので, もっと単純に \[ x = -\frac{1}{\lambda}\ln u\quad (u \sim U(0,1)) \] によって生成すれば十分です.

以下がPythonでの実装例です.

from math import log
from random import random
def exponential(lam):
    return -log(random())/lam

\(\lambda = 5\) として前頁の関数で生成した乱数サンプルを 1000 個作ってヒストグラム化したもの. (幅0.1の区間毎に頻度を計算.)

応用:ポアソン分布

指数分布に従う乱数を生成する事が出来たので, ポアソン分布に従う乱数も生成する事が出来ます.

復習すると, ある事象の単位時間あたりの発生回数がポアソン分布, その発生間隔が指数分布に従うのでした.

従って, 一様乱数 \(u_1,u_2,\ldots,u_n\) を次々に生成して \[ \begin{aligned} &-\frac{1}{\lambda}\ln u_1 -\frac{1}{\lambda}\ln u_2 -\cdots -\frac{1}{\lambda}\ln u_n \leq 1 \\ \Leftrightarrow & \color{yellow}{ e^\lambda u_1u_2\cdots u_n \geq 1 } \end{aligned} \] を満たす, 最大の \(n\) を返せば良いです.

以下が実装例です.

from random import random
from math import exp
def poisson(lam):
    v = exp(lam) * random()
    k = 0
    while v >= 1:
        v *= random()
        k += 1
    return k

逆変換法(離散的分布)

離散的な分布についても同様の事が出来ます.

\(U\) が標準一様分布に従うとして, 求める確率変数 \(X\) を \[ X = \left\{\begin{array}{cl} 1 & (0\leq U < 1/6) \\ 2 & (1/6\leq U < 2/6) \\ 3 & (2/6\leq U < 3/6) \\ 4 & (3/6\leq U < 4/6) \\ 5 & (4/6\leq U < 5/6) \\ 6 & (5/6\leq U < 1) \end{array}\right. \] と定めれば良いです.

以下がPythonでのコーディング例です.

from random import random

def dice():
    r = random()
    for i in range(1,7):
        if r < i/6.0: return i

【注】
実際には擬似乱数は離散的な値しか取らない為, この方法では均等に\(1/6\)の確率で出現するサイコロにはなりません. 完全に均等にしたい場合には, 以下の様な手段を取ります.

• 整数値の乱数を振る.(仮に32bitとすると0以上4294967295以下の値を得る.)
• 値域が6の倍数個となるように, 4294967291以上の値は捨てて振り直す.
• (6で割った余り)+1 を返す.

棄却法

逆関数の計算が困難な場合には, 以下の方法があります.

乱数の組 \((u_1,u_2)\) は下図の長方形内の１点に対応します. 棄却法では, \(y=\rho(x)\) で囲まれた領域に点が入った場合だけ採用するという方法を取ります.

\[ \rho'(x) = 30(1-x)^4 - 120x(1-x)^3 = 30(1-5x)(1-x)^3 \] なので, \(\rho(x)\) は \(x = \frac{1}{5}\) で最大値約\(2.5\) を取ります. よって \(c = 2.6\) と置けば十分です.

from random import random
def beta25():
    while True:
        u1 = random()
        u2 = random() * 2.6
        if u2 <= 30*u1*(1-u1)**4:
            return u1

1000サンプルの乱数を生成して頻度をプロットしたのが以下の図です.

先ほどの単純な棄却法では, \(y=\rho(x)\)の外の部分(棄却域)の面積が大きいときに反復回数が増えてしまいます.

これは, 既知の確率密度関数\(y=\varphi(x)\)を上手く選んで, それを適当に縦の伸ばしたものを矩形の代わりに利用する事で改善できます.

正規分布の乱数生成法ですが, これは必要な知識をまだ説明していないので方法だけ紹介します.

推定

本節では, 調べたい集団についてそのサンプルに基いて推定を行う方法について解説します.

調べる対象の集団(確率分布)を母集団, 母集団からランダムに取り出したデータの集合を標本と呼びます.

母集団の特徴量は 母平均,母分散,母標準偏差 などの様に呼びます. 一般に,これらは母数とも呼ばれます. また, 標本の特徴量は 標本平均,標本分散,標本標準偏差 などの様に呼びます.

点推定

「母平均は \(\mu\) である」といった形で, その値をピンポイントで推定する事を 点推定 と言います.

点推定には以下の方法があります.

• 不偏推定量による推定
• 最尤推定量による推定

不偏推定量による推定

【証明】
母平均を \(\mu\) とすると 各標本について \(E[X_i] = \mu\) であるから \[ \begin{aligned} &E\left[\frac{X_1+X_2+\cdots +X_n}{n}\right] = \frac{1}{n}(E[X_1]+E[X_2]+\cdots + E[X_n]) \\ &= \frac{\mu+\mu+\cdots+\mu}{n} = \mu \end{aligned} \] □

\( (\text{標本数})-1 \) で割り算をしている事に注意してください.

【証明】
母平均を \(mu\), 母分散を \(\sigma^2\) とおくと, \(E[X_i]=\mu, V[X_i]=\sigma^2\) である. また \[ s^2 = \frac{1}{n-1}\sum_{i=1}^n(X_i-\overline{X})^2 \] とおく. ここで \[ \begin{aligned} (X_i - \overline{X})^2 &= \{(X_i-\mu)-(\overline{X}-\mu)\}^2 \\ &(X_i-\mu)^2 + (\overline{X}-\mu)^2 -2(X_i-\mu)(\overline{X}-\mu) \end{aligned} \] より \[ \begin{aligned} \sum_{i=1}^n(X_i-\overline{X})^2 &= \sum_{i=1}^n (X_i-\mu)^2 + n(\overline{X}-\mu)^2 - 2n(\overline{X}-\mu)^2 \\ &= \sum_{i=1}^n (X_i-\mu)^2 - n(\overline{X}-\mu)^2 \end{aligned} \] である.

ここで, \(E[\overline{X}]=\mu\) である事に注意すれば \[ \begin{aligned} E[(\overline{X}-\mu)^2] &= V[\overline{X}] = V\left[\frac{1}{n}(X_1+X_2+\cdots+X_n)\right] \\ &= \frac{1}{n^2}(V[X_1]+V[X_2]+\cdots+V[X_n]) = \frac{1}{n}\sigma^2 \end{aligned} \] となるから, \[ \begin{aligned} E[s^2] &= \frac{1}{n-1}\left\{ \sum_{i=1}^n E[(X_i-\mu)^2] - n\cdot \frac{\sigma^2}{n} \right\} \\ &= \frac{1}{n-1}(n\sigma^2 - \sigma^2) \\ &= \sigma^2 \end{aligned} \] □

実験をしてみましょう. パラメータ \(\lambda\) のポアソン分布(母平均も母分散も\(\lambda\))に従う乱数を100個生成して母平均・母分散を推定してみます.

>>> import numpy as np
>>> N = 1000
>>> data = np.random.poisson(5, N)
>>> np.average(data)
4.9359999999999999
>>> np.var(data, ddof=1)  # 不偏分散
4.7466506506506221
>>> np.var(data)          # 標本分散
4.7419039999999715

最尤推定量による推定

最尤推定量とは, 標本の得られる確率が最大となる推定量です.

18回1の目が出る確率は \[ P = \binom{100}{18}p^{18}(1-p)^{82} \] です. これを \(p\) で微分すれば \[ \begin{aligned} \frac{\mathrm{d}P}{\mathrm{d} p} &= \binom{100}{18}\left\{18p^17(1-p)^{82}-82p^{18}(1-p)^{81}\right\} \\ &= \binom{100}{18}p^{17}(1-p)^{81}(18-100p) \end{aligned} \] となるので, \(p = \frac{18}{100} = 0.18\)が最尤推定量となります.

同時確率関数・同時密度関数は積の形になる事が多いので, 対数を取って利用する事が多いです. 対数関数は単調増加するので, 対数尤度が最大となる \(\theta\) は最尤推定量と同じです.

【証明】
母平均を\(\mu\),母分散を\(\sigma^2\) とすると, 同時確率密度関数は \[ \rho(x_1,x_2,\ldots,x_n) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{-\frac{(x_i-\mu)}{2\sigma^2}\right\} \] であるから, 対数尤度関数は \[ \begin{aligned} \ln L(\mu,\sigma^2) &= \sum_{i=1}^n \ln\left[\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{-\frac{(x_i-\mu)}{2\sigma^2}\right\} \right] \\ &= n\ln\frac{1}{\sqrt{2\pi\sigma^2}} - \frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-\mu)^2 \end{aligned} \] となる. これを \(\mu,\sigma^2\) で偏微分すれば \[ \begin{aligned} \frac{\partial }{\partial \mu} \ln L(\mu,\sigma^2) &= \frac{1}{\sigma^2}\sum_{i=1}^n (x_i-\mu) \\ \frac{\partial}{\partial(\sigma^2)}\ln L(\mu,\sigma^2) &= -\frac{n}{2\sigma^2}+\frac{1}{2\sigma^4}\sum_{i=1}^n(x_i-\mu)^2 \end{aligned} \]

従って, 対数尤度が最大となるのは \[ \begin{aligned} \mu &= \frac{1}{n}\sum_{i=1}^n x_i \\ \sigma^2 &= \frac{1}{n}\sum_{i=1}^n(x_i-\mu)^2 \end{aligned} \] のときである. □

【証明】
対数尤度関数は \[ \ln L(\lambda) = \ln\prod_{i=1}^n \lambda e^{-\lambda x_i} = \ln\left(\lambda^n e^{-\lambda \sum_{i=1}^n x_i}\right) = n\ln\lambda-\lambda\sum_{i=1}^nx_i \] であるので, \[ \frac{\mathrm{d}}{\mathrm{d}\lambda}\ln L(\lambda) = \frac{n}{\lambda} - \sum_{i=1}^n x_i \] となる. よって, \(L(\lambda)\) が最大となるのは, \[ \frac{n}{\lambda}-\sum_{i=1}^n x_i=0\Leftrightarrow \lambda=\frac{n}{\sum_{i=1}^nx_i} \] のときである. □

標本平均は\( \frac{0.5+34.7+5.0+13.5+11.7+7.0+0.7+8.1+26.7+9.0}{10}=11.69\)なので, \(\lambda\) の最尤推定量は\(1/11.69\approx 0.086\) となります. 指数分布では, 1単位時間あたり平均して\(\lambda\)回イベントが発生するのですから, 1秒あたりの平均アクセス数は約86回と推定する事が出来ます.

【証明】
対数尤度関数は \[ \begin{aligned} \ln L(\lambda) &= \ln\left(\frac{\lambda^{x_1}e^{-\lambda}}{x_1!} \frac{\lambda^{x_2}e^{-\lambda}}{x_2!} \cdots \frac{\lambda^{x_n}e^{-\lambda}}{x_n!}\right)\\ &= -\ln(x_1!x_2!\cdots x_n!)+(\sum_{i=1}^nx_i)\ln\lambda-n\lambda \end{aligned} \] であるので, \[ \frac{\mathrm{d}}{\mathrm{d}\lambda}\ln L(\lambda) = \frac{\sum_{i=1}^n x_i}{\lambda} - n \] となる. よって, \( L(\lambda) \) が最大となるのは, \[ \frac{\sum_{i=1}^n x_i}{\lambda} - n =0\Leftrightarrow \lambda=\frac{\sum_{i=1}^n x_i}{n} \] の時である. □

区間推定

「母平均は95%の確率で～から～の範囲にある」といった形で, 特徴量の値の取る範囲とその確率を推定する事を 区間推定 と言います.

\(\alpha\)の値としては例えば \(\alpha=0.05\) などがよく利用されます. この場合, 母数 \(\theta\) の95%信頼区間とは \[ P(a\leq \theta \leq b) = 0.95 \] となる \([a,b]\) という事になります.

標準正規分布は左右対称なので, 上側 \(\alpha\) 点の値が計算出来れば両側 \(\alpha\) 点も計算出来ます. よく使いそうな値だけ下に書いておきます. \[ \begin{array}{|c|c|} \hline \alpha & \text{上側$\alpha$点} \\ \hline 0.1 & 1.281552 \\ 0.05 & 1.644854 \\ 0.025 & 1.959964 \\ \hline \end{array} \]

使っている言語でのこれらの計算方法を確認しておきましょう. 以下は SciPyでの上側 \(\alpha\) 点,両側\(\alpha\)点の計算例です.

>>> from scipy import stats
>>> stats.norm.ppf(0.9)   # 上側10%点 = 両側20%点
1.2815515655446004
>>> stats.norm.ppf(0.95)  # 上側5%点 = 両側10%点
1.6448536269514722

母集団が\(N(\mu,\sigma^2\)\)に従っているならば, 標本平均 \[ \overline{X}=\frac{X_1+X_2+\cdots+X_n}{n} \] は\( \sim N(\mu,\sigma^2/n) \)に従います. よって \[ Y = \frac{\sqrt{n}}{\sigma}(\overline{X}-\mu) \] は標準正規分布 \(N(0,1)\) に従いますので, 両側\(\alpha\)点を\(z_\alpha\)とすれば \[ P(|Y|\leq z_\alpha)=1-\alpha\Leftrightarrow P\left(|\overline{X}-\mu|\leq \frac{\sigma}{\sqrt{n}}z_\alpha\right)=1-\alpha \] となります.

\(\alpha=0.05\)の両側\(\alpha\)点の値は \(z_\alpha \approx 1.959964\) であるので 95% 信頼区間は \[ 3 - \frac{\sqrt{10}}{\sqrt{100}}1.959964 \leq \mu \leq 3 + \frac{\sqrt{10}}{\sqrt{100}}1.959964 \] つまり, 有効数字を3桁とすると\(3-0.62 \leq \mu \leq 3+0.62\)となります.

続いて,母分散が未知の場合を説明します. \(t\)分布という分布が必要になります.

\(t\) 分布の確率密度関数の具体的な表式の導出にはガンマ関数などの知識が必要となるため省略します. 以下が自由度\(1,2,3\) の\(t\)分布の密度関数のグラフとなります.

\(t\)分布に対しても上側\(\alpha\)点,両側\(\alpha\)点が定められます. \(t\)分布も左右対称の密度関数を持つので, 上側\(\alpha\)点が計算出来れば十分です.

以下は 上側 \(\alpha\) 点の表です.

既に述べましたが, 中心極限定理によって正規分布以外の分布に対しても信頼区間の計算を行う事が出来ます.

例えば, 以下のような問題を考えてみましょう.

母集団がベルヌーイ分布である場合の信頼水準\(1-\alpha\)の信頼区間を計算すれば良いです.
ベルヌーイ分布のパラメータを\(p\)とすれば, 母平均は\(p\),母分散は\(p(1-p)\)ですから, 母分散が既知の場合の結果を使う事が出来ます. \(n\) が十分大きいとき, 標本平均を\(\overline{X}\), 正規分布の両側\(\alpha\)点を\(z_\alpha\) とすると 信頼水準\(1-\alpha\)の信頼区間は \[ \overline{X} - \sqrt{\frac{p(1-p)}{n}}z_\alpha \leq p \leq \overline{X} + \sqrt{\frac{p(1-p)}{n}}z_\alpha \] となります. この両辺には未知数\(p\)が含まれてしまっているので, \(p\)について解き直します.

\[ \left| p-\overline{X} \right| ^2 \leq \frac{p(1-p)}{n}z_\alpha^2 \] と同値変形して, 両辺を \(p\) について整理すれば \[ \left(1+\frac{z_\alpha^2}{n}\right)p^2 - 2\left(\overline{X}+\frac{z_\alpha^2}{2n}\right)p + \overline{X}^2 \leq 0\] となります.

これを解いて \[ \frac{\overline{X}+\frac{z_\alpha^2}{2n}-z_\alpha\sqrt{\frac{\overline{X}(1-\overline{X})}{n}+\frac{z_\alpha^2}{4n^2}}}{1+\frac{z_\alpha^2}{n}} \leq p \leq \frac{\overline{X}+\frac{z_\alpha^2}{2n}+z_\alpha\sqrt{\frac{\overline{X}(1-\overline{X})}{n}+\frac{z_\alpha^2}{4n^2}}}{1+\frac{z_\alpha^2}{n}} \] が求める信頼区間となります. \(n\)が十分大きい場合には \(\frac{z_\alpha^2}{n}\) を無視する事ができるので, \[ \color{yellow}{ \overline{X} - \sqrt{\frac{\overline{X}(1-\overline{X})}{n}}z_\alpha \leq p \leq \overline{X} \sqrt{\frac{\overline{X}(1-\overline{X})}{n}}z_\alpha} \] となります.

先ほどの問題に適用してみましょう.

以下の計算より, 求める95%信頼区間は \[ 0.2 - 0.03 \leq p \leq 0.2 + 0.03 \] つまり \(\color{yellow}{0.17\leq p \leq 0.23}\) となります.

>>> from scipy import stats
>>> from math import sqrt
>>> x_bar = 120.0 / 600
>>> z_alpha = stats.norm.ppf(0.975)
>>> sqrt(x_bar*(1-x_bar)/600)*z_alpha
0.032006077842368738
>>> x_bar
0.2

母分散の区間推定の為には以下の分布が必要になります.

\(\chi^2\)分布の具体的な立式は省略します. 自由度 \(1,2,3\) についての密度関数のグラフは以下のようになっています. 正規分布・\(t\)分布と違って左右対称ではないので注意してください.

以下の計算より \(5.2 \leq \sigma^2 \leq 9.1\) となります.

>>> from scipy import stats
>>> left = stats.chi2.ppf(0.975, 99)    # 上側α/2点
>>> right = stats.chi2.ppf(0.025, 99)   # 上側1-α/2点
>>> 100*6.7 / left
5.2171750887470614
>>> 100*6.7 / right
9.1329080522400083

今回はここで終わります。

次回は検定, 標本数の決定について説明します. 時間があれば多変量分布についても解説します.