2017-06-26

二項分布の極限が正規分布になることの証明【正規近似の証明】

数学統計

二項分布の試行回数を無限大に大きくしていくと正規分布に近づくことが知られている。
しかし、その証明は意外と知られていない。（中心極限定理でも証明は可能ではあるが、回りくどすぎて初学者には不親切）

$e$ や $\pi$ が出てくるのがどうしても不思議で、自分なりに解こうとしたり調べてみた。
ようやくわかってきたので、自分的にわかりやすかった方法を記録しておく。
ただ、これはわかりやすさ重視の証明法らしい。
数学的に厳密じゃないところとか突っ込んでくれると嬉しいです。

予備知識

二項分布 : $\displaystyle B(n,p)={}_n C_{x} p^x(1-p)^{n-x}$

正規分布 : $\displaystyle \mathcal{N}(\mu,\,\sigma^{2}) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(\mu - x)^2}{2\sigma^2}\right)$

二項分布の正規近似 : $\displaystyle \lim_{n \to \infty} B(n,p)=\mathcal{N}(np,\,np(1-p))$

証明

方針

二項分布の対数をとった関数をTaylor展開し二次の項まで近似すると正規分布の対数になることを示す。

証明

二項分布の自然対数をとる．
$\displaystyle \begin{equation} \begin{split} \log{B(n,p)} &= \log{{}_n C_{x} p^x(1-p)^{n-x}}\\ &= \log{\frac{n!}{x!(n-x)!} p^x(1-p)^{n-x}}\\ &= \log{n!} - \log{x!} - \log{(n-x)!} + x\log{p} + (n-x)\log{(1-p)}　 \end{split} \end{equation}$
ここでTaylor展開を考えるにあたり，
$\log{x!},\, \log{(n-x)!}$
の項について考える．

補題として，
[補題１]
充分に大きい自然数 $k$ に対して
$\displaystyle \log{k!} \simeq \int_1^k\log{t} dt$
が成り立つことを示す。
[補題１証明]
2以上の自然数 $a$ に対して， $\log{(t)}$ は単調増加であるので
$\displaystyle \log{(a-1)} < \int_{a-1}^a \log{t} dt < \log{a}$
が成り立つ．
したがって
$\displaystyle \begin{equation} \log{(k-1)!} < \int_1^{k} \log{t} dt < \log{k!} \\ \frac{\log{(k-1)!}}{\log{k!}} < \frac{\int_1^{k} \log{t} dt}{\log{k!}} < 1 \end{equation}$

十分に大きい $k$ に対して
$\displaystyle \frac{\log{(k-1)!}}{\log{k!}} = 1 - \frac{\log{k}}{\log{k!}} \simeq 1$
である．*1
上式からはさみうちの原理より、
$\displaystyle \frac{\int_1^{k} \log{t} dt}{\log{k!}} \simeq 1$
となる．
よって
$\displaystyle \log{k!}　\simeq \int_1^{k} \log{t} dt$
[補題１証明終]

補題１の結果を元の式に代入数する．
$x ,\, (n-x)$ が十分に大きいとき
$\displaystyle \begin{equation} \begin{split} \log{B(n,p)} &= \log{n!} - \log{x!} - \log{(n-x)!} + x\log{p} + (n-x)\log{(1-p)}\\ &\simeq \log{n!} - \int_1^{x} \log{t} dt - \int_1^{n-x} \log{t} dt + x\log{p} + (n-x)\log{(1-p)} \end{split} \end{equation}$
となる．

この変形によって微分を考えることができて，
一階微分、二階微分はそれぞれ
$\displaystyle \begin{equation} \begin{split} \\(\log{B(n,p)})' &\simeq - \log{x} + \log{(n-x)} + \log{p} - \log{(1-p)} \\ \\(\log{B(n,p)})'' &\simeq - \frac{1}{x} - \frac{1}{n-x} \end{split} \end{equation}$
となる．

$\log{B(n,p)}$ を平均値の周りでTaylor展開する．
ただし、平均値と分散は二項分布の性質から定まる（平均 $\mu = np$ 、分散 $\sigma^2 = np(1-p)$ ）

$\displaystyle \begin{equation} \begin{split} \log{B(n,p)} &\simeq \log{B(x=\mu | n,p)} + (\log{B(x=\mu | n,p)})'(x-\mu) + \frac{1}{2!}(\log{B(x=\mu | n,p)})''(x - \mu)^2 \\ &\simeq \log{B(x=\mu | n,p)} - \frac{1}{2np(1-p)}(x - \mu)^2 \\ &= \log{B(x=\mu | n,p)} - \frac{(x - \mu)^2}{2\sigma^2} \end{split} \end{equation}$

ここで，第一項は定数なので、 $\log{C}$ とおいておく．
（近似の過程で確率密度の制約条件「定義域全域で積分すると１になる」が抜けてしまうため、改めて係数を求める必要がある．）

求める確率分布を $N(x)$ とすると，
$\displaystyle \begin{equation} \begin{split} \log{N(x)} &\simeq \log{C} - \frac{1}{2\sigma^2}(x - \mu)^2 \\ N(x) &\simeq C \exp{\left(-\frac{(x - \mu)^2}{2\sigma^2}\right)} \end{split} \end{equation}$
これで変数部分まで求めることができた．

最後に，
先ほども述べたように，確率密度の制約条件を満たすよう係数を決定しなければいけない．
すなわち，
$\displaystyle \int_{-\infty}^{\infty} N(x) dx = 1$
となる $C$ を求める．

補題として、
[補題2]
$\displaystyle \int_{0}^{\infty} e^{-x^2} dx = \frac{\sqrt{\pi}}{2}$
を示す．（※ガウス積分）
[補題2証明]
$\displaystyle I(R) = \int_{0}^{R} e^{-x^2} dx$
とおくと、証明すべきは
$\displaystyle \lim_{R \to \infty}I(R) = \frac{\sqrt{\pi}}{2}$
である．

$\displaystyle (I(R))^2 = I(R)I(R) = \int_{0}^{R}\int_{0}^{R} e^{-(x^2+y^2)}dxdy$
とおくと，この積分領域は $xy$ 平面上における一辺 $R$ の正方形である．
一方で次のような扇状の積分範囲を定義すると，
$D(R) = \left\{(x,y)\,|\, x\geq 0,\, y \geq 0,\, x^2+y^2 \le R^2 \right\}$
次のような大小関係が成り立つ
$\displaystyle \int\int_{D(R)} e^{-(x^2+y^2)}dxdy < (I(R))^2 < \int\int_{D(\sqrt{2}R)} e^{-(x^2+y^2)}dxdy$

ここで，
$\displaystyle J(R) = \int\int_{D(R)} e^{-(x^2+y^2)}dxdy$
とする．
$\displaystyle \lim_{R \to \infty} J(R) = \frac{\pi}{4}$
となることから，はさみうちの原理より $(I(R))^2 = \sqrt{\pi}/4$ となることを示す．

$J(R)$ を極座標表示に変換して計算する
$x = r\cos{\theta} , \, y = r\sin{\theta} \\ 0 \leq r \leq R ,\, 0 \leq \theta \leq \frac{\pi}{2}$
として，
$\displaystyle \begin{equation} \begin{split} J(R) &= \int_{0}^{R}\int_{0}^{\pi/2} e^{-r^2} r d\theta dr \\ &= \left\{\int_{0}^{R} re^{-r^2} dr \right\} \left\{\int_{0}^{\pi/2} 1 d\theta \right\} \\ &= \frac{\pi}{2}\left[ - \frac{1}{2} e^{-r^2} \right]_{0}^{R} \\ &= \frac{\pi}{4}\left( 1 - e^{-R^2} \right) \end{split} \end{equation}$

$R$ が十分に大きいとき，
$\displaystyle \begin{equation} J(R) = \frac{\pi}{4} \\ J(\sqrt{2}R) = \frac{\pi}{4} \end{equation}$
であることがわかる．

はさみうちの原理より，
$\displaystyle \begin{equation} \begin{split} (I(R))^2 &= \frac{\pi}{4} \\ I(R) &= \frac{\sqrt{\pi}}{2} \end{split} \end{equation}$

が示された．
$\displaystyle \left(\int_{0}^{\infty} e^{-x^2} dx = \frac{\sqrt{\pi}}{2} \right)$
[補題2証明終]

補題２の結果から，係数 $C$ を求める．
$x = \sqrt{2}\sigma X +\mu$
と変数変換することで，
$\displaystyle \begin{equation} \begin{split} \int_{-\infty}^{\infty} N(x) dx &= \int_{-\infty}^{\infty} C \exp{\left(-\frac{(x - \mu)^2}{2\sigma^2}\right)} dx \\ &= \sqrt{2}\sigma C \int_{-\infty}^{\infty} e^{-X^2}dX \\ &= \sqrt{2\pi\sigma^2} C \end{split} \end{equation}$
から，
$\displaystyle C=\frac{1}{\sqrt{2\pi\sigma^2}}$
が求まる．

以上から，
$\displaystyle \lim_{n \to \infty} B(n,p) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$
となり，二項分布の極限は正規分布に近似されることが証明された．
[証明終]

*1:コメントに証明のあらすじを書いてくれている方がいます。ありがとうございます。

2017-06-23

直下のファイルを走査し一覧表示する [bash, python]

あるディレクトリの直下のファイルに対して巡回しながら処理したい時
シェルスクリプトとpythonでどうするか。

bash

files="./hoge/*"
for filepath in $files; do
  echo $filepath
done

./hoge/1.txt
./hoge/2.txt
のようにディレクトリごと表示される
ファイル名だけ表示したい場合はこう書けば良い

${filepath##*/}

補足

zipファイルの中にzipファイルがあるもの、例えば

a.zip
- nbghuio.zip
- bghjkoi.zip
b.zip
- bhujiko.zip
- qvdshjt.zip

こんな感じのものをまとめて解凍したい時

#!/bin/bash
for zipfile `\find . -name '*.zip'` ; do
    unzip -q ${zipfile}
done

for i in "./*" ; do
    for file in `\find $i -name '*.zip'`; do
        unzip -q -o ${file}
    done
done

こんなスクリプトになる。

Python

import os

directory = "./hoge/"
for root, dirs, files in os.walk(directory):
    for file in files:
        print(file)
        print(root+file)

fileだけの場合は
1.txt
のようになる。rootにディレクトリの文字列があるので
./hoge/1.txt
のようにしたい場合はroot+fileのように書く

2017-06-22

pythonで完全数と友愛数と婚約数を求める

数学

以前の記事の続きとして、せっかくなので友愛数・完全数・婚約数を求めてみたい。

約数関数（再掲）

divisor(n)として約数関数を使うので、再掲しておく。
sympyを使う。

import sympy
def divisor(n):
    factors = sympy.factorint(n)
    rst = 1
    for i,j in factors.iteritems():
        rst *= (pow(i,j+1) - 1)/(i-1)
    return rst

完全数

完全数（かんぜんすう，英: perfect number）とは、自分自身を除く正の約数の和に等しくなる自然数のことである。完全数の最初の3個は 6 (= 1 + 2 + 3)、28 (= 1 + 2 + 4 + 7 + 14)、496 (= 1 + 2 + 4 + 8 + 16 + 31 + 62 + 124 + 248) である。「完全数」は「万物は数なり」と考えたピタゴラスが名付けた数の一つであることに由来する[1]が、彼がなぜ「完全」と考えたのかについては何も書き残されていないようである
wikipedia:完全数

max = 100000
rst = []
for n in range(1,max+1):
    if divisor(n) == 2*n:
        rst.append(n)

結果

rst = 
[6, 28, 496, 8128]

友愛数

友愛数（ゆうあいすう、英: amicable numbers）とは、異なる 2 つの自然数の組で、自分自身を除いた約数の和が、互いに他方と等しくなるような数をいう。親和数（しんわすう）とも呼ばれる。
最小の友愛数の組は (220, 284) である。
220 の自分自身を除いた約数は、1, 2, 4, 5, 10, 11, 20, 22, 44, 55, 110 で、和は 284 となる。一方、284 の自分自身を除いた約数は、1, 2, 4, 71, 142 で、和は 220 である。
wikipedia:友愛数

max = 100000
rst = []
for n in range(1,max+1):
    pair = divisor(n) - n
    if n < pair and pair <= max:
        if divisor(pair) - pair == n:
            rst.append([n,pair])

結果

rst = 
[[220, 284],
 [1184, 1210],
 [2620, 2924],
 [5020, 5564],
 [6232, 6368],
 [10744, 10856],
 [12285, 14595],
 [17296, 18416],
 [63020, 76084],
 [66928, 66992],
 [67095, 71145],
 [69615, 87633],
 [79750, 88730]]

婚約数

婚約数（こんやくすう、英: betrothed numbers）とは、異なる 2 つの自然数の組で、1 と自分自身を除いた約数の和が、互いに他方と等しくなるような数をいう。準友愛数（じゅんゆうあいすう、quasi-amicable numbers）とも呼ばれる。
最小の婚約数の組は (48, 75) である。
48 の 1 と自分自身を除いた約数は、2, 3, 4, 6, 8, 12, 16, 24 で、和は 75 となる。一方、75 の 1 と自分自身を除いた約数は、3, 5, 15, 25 で、和は 48 である。
wikipedia:婚約数

max = 100000
rst = []
for n in range(1,max+1):
    pair = divisor(n) - n - 1
    if n < pair and pair <= max:
        if divisor(pair) - pair - 1 == n:
            rst.append([n,pair])

結果

rst = 
[[48, 75],
 [140, 195],
 [1050, 1925],
 [1575, 1648],
 [2024, 2295],
 [5775, 6128],
 [8892, 16587],
 [9504, 20735],
 [62744, 75495]]

技術メモ

役に立てる技術的な何か、時々自分用の覚書。幅広く色々なことに興味があります。