連休初日から2日連続有給ぶちかましたkimotonですこんにちは。

中心極限定理とは

大数の法則によると、ある母集団から無作為抽出した標本の平均は標本の大きさを大きくすると母平均に近づく。これに対し中心極限定理は標本平均と母平均との誤差を論ずるものである。多くの場合、母集団の分布がどんな分布であっても、その誤差は標本の大きさを大きくしたとき近似的に正規分布に従う。

母集団が どんな 分布であっても、その母集団の平均を $\mu$ 、分散を $\sigma ^2$ とした時、そこから $n$ 個ランダムサンプリングした標本の平均は、 $n$ が十分に大きい時、近似的に正規分布に従うことがわかっているのです。

今回はこの中心極限定理について、n=1の標本を用いて簡単に実証してみたいと思います。

正規分布とは

ガウス分布、誤差分布とも呼ばれ、平均値の付近に集積するようなデータ分布を表した連続的な変数に関する確率分布です。確率密度関数は下記の通りです。

${{\displaystyle f(x)={\frac {1}{\sqrt {2\pi \sigma^{2}}}}\exp !\left(-{\frac {(x-\mu )^{2}}{2\sigma^{2}}}\right)\quad (x\in \mathbb {R} )}}$

pythonで実装する際はこんな感じ。そのままです。

def normal_pdf(x, mu=0, sigma=1):
    sqrt_two_pi = math.sqrt(2 * math.pi)
    return (math.exp(-(x-mu) ** 2 / 2 / sigma ** 2) / (sqrt_two_pi * sigma))

$(\mu, \sigma) = (0, 1), (0, 2), (0, 0.5), (-1, 1)$ について正規分布を描画

import matplotlib.pyplot as plt
xs = [x / 10.0 for x in range(-50, 50)]
plt.plot(xs,[normal_pdf(x,sigma=1) for x in xs],'-',label='mu=0,sigma=1')
plt.plot(xs,[normal_pdf(x,sigma=2) for x in xs],'--',label='mu=0,sigma=2')
plt.plot(xs,[normal_pdf(x,sigma=0.5) for x in xs],':',label='mu=0,sigma=0.5')
plt.plot(xs,[normal_pdf(x,mu=-1)   for x in xs],'-.',label='mu=-1,sigma=1')
plt.legend()
plt.title("Various Normal pdfs")
# plt.show()

f:id:kimoppy126:20190716125312p:plain

正規分布の累積分布関数

累積分布関数の差を求めることで、離散確率分布の確率を求めることができます。
正規分布の累積分布関数は誤差関数は深い繋がりがあります。

誤差関数は下記
$\displaystyle \mathrm{erf}(x):=\frac{2}{\sqrt{\pi}} \int_{0}^{x} \exp(-t) dt$

正規分布の累積分布関数はこれを使って下記のように表せる。
$\displaystyle F(x) = \frac{1}{2} \{ 1 + \mathrm {erf}(\frac{x-\mu}{\sqrt{2}\sigma}) \}$

pythonで実装する際はこんな感じ。

def normal_cdf(x, mu=0,sigma=1):
    return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2

$(\mu, \sigma) = (0, 1), (0, 2), (0, 0.5), (-1, 1)$ について、正規分布の累積密度分布を描画

xs = [x / 10.0 for x in range(-50, 50)]
plt.plot(xs,[normal_cdf(x,sigma=1) for x in xs],'-',label='mu=0,sigma=1')
plt.plot(xs,[normal_cdf(x,sigma=2) for x in xs],'--',label='mu=0,sigma=2')
plt.plot(xs,[normal_cdf(x,sigma=0.5) for x in xs],':',label='mu=0,sigma=0.5')
plt.plot(xs,[normal_cdf(x,mu=-1) for x in xs],'-.',label='mu=-1,sigma=1')
plt.legend(loc=4) # bottom right
plt.title("Various Normal cdfs")
plt.show()

f:id:kimoppy126:20190716194207p:plain

二項分布で中心極限定理

二項分布は、ベルヌーイ試行を独立に行ったときの成功回数を確率変数とする離散確率分布です。
成功確率 $p$ 、試行回数 $N$ としたとき、二項分布は下記の確率密度関数で表されます。

$\displaystyle Bin(x|N, p) = {}_N C_x p^{x}(1-p)^{N-x}$

これをpythonで実装する際はこんな感じ。

def binomial_pdf(p, n, r):
    return (p ** r) * ((1-p)**(n-r)) * math.factorial(n) / (math.factorial(n - r) * math.factorial(r))

$(n, p) = (100, 0.25), (100, 0.5), (100, 0.75)$ の二項分布を描画

import matplotlib.pyplot as plt

xs = [x for x in range(0, 100)]
plt.plot(xs,[binomial_pdf(n=100, p=0.25, r=x) for x in xs],'--',label='p=0.25,n=100')
plt.plot(xs,[binomial_pdf(n=100, p=0.5, r=x) for x in xs],':',label='p=0.5,n=100')
plt.plot(xs,[binomial_pdf(n=100, p=0.75, r=x) for x in xs],'-.',label='p=0.75,n=100')
plt.legend()
plt.title("Various Binomial pdfs")
plt.show()

f:id:kimoppy126:20190716230020p:plain

二項分布からの無作為抽出はベルヌーイ試行n回の合計を求めればよいため、下記の通り定義できます。

# ベルヌーイ試行
def bernoulli_trial(p):
    return 1 if random.random() < p else 0
# 二項分布
def binomial(n, p):
    return sum(bernoulli_trial(p) for _ in range(n))

二項分布を母集団とするn=10000の標本が正規分布に従うことをヒストグラムを用いて可視化することで確かめます。

def make_hist(p, n, num_points):
    data = [binomial(n, p) for _ in range(num_points)]
    histogram = Counter(data)
    plt.bar([x - 0.4 for x in histogram.keys()],
    [v / num_points for v in histogram.values()],
    0.8,
    color='0.75')
    mu = p * n
    sigma = math.sqrt(n * p * (1 - p))

    # use a line chart to show the normal approximation
    xs = range(min(data), max(data) + 1)
    ys = [normal_cdf(i + 0.5, mu, sigma) - normal_cdf(i - 0.5, mu, sigma) for i in xs]
    plt.plot(xs,ys)
    plt.title("Binomial Distribution vs. Normal Approximation")
    plt.show()

f:id:kimoppy126:20190716181422p:plain

その他

numpy.randomパッケージでは、様々な分布からランダムサンプリングを行う関数を使用することができます。
以下では、この関数を使用して各分布からの無作為抽出値が正規分布に近似することを確かめます。

正規分布で中心極限定理

from numpy.random import normal
from numpy import arange
def make_normal_hist(mu, sigma, num_points):
    data = [normal(mu, sigma) for _ in range(num_points)]
    plt.hist(data,
    bins=30,
    rwidth = 0.8,
    color='0.75',
    normed=True)

    # use a line chart to show the normal approximation
    xs = arange(min(data), max(data), 0.01)
    ys = [normal_pdf(i, mu, sigma) for i in xs]
    plt.plot(xs,ys)
    plt.title("Poisson Distribution vs. Normal Approximation")
    plt.show()
    plt.title("Distribution(n={}) from Normal Distribution".format(num_points))

n=100

> make_normal_hist(100, 1, 100)

f:id:kimoppy126:20190718073203p:plain

n=10000

> make_normal_hist(100, 1, 10000)

f:id:kimoppy126:20190718073135p:plain

n数が増えるに従いより正規分布に近い分布を示すことがわかります。

ポアソン分布で中心極限定理

from numpy.random import poisson
def make_poisson_hist(lam, num_points):
    data = [poisson(lam) for _ in range(num_points)]
    histogram = Counter(data)
    plt.bar([x - 0.4 for x in histogram.keys()],
    [v / num_points for v in histogram.values()],
    0.8,
    color='0.75')
    mu = lam
    sigma = math.sqrt(lam)

    # use a line chart to show the normal approximation
    xs = range(min(data), max(data) + 1)
    ys = [normal_cdf(i + 0.5, mu, sigma) - normal_cdf(i - 0.5, mu, sigma) for i in xs]
    plt.plot(xs,ys)
    plt.title("Distribution(n={}) from Poisson Distribution".format(num_points))
    plt.show()