定理

$n$ 個の確率変数 $X_1, \cdots, X_n$ が

• 互いに独立している
• 平均 $\mu$、分散 $\sigma^2$ の同一の確率分布（正規分布でなくて良い）に従う

とき、これらの平均

$$\bar{X} = X_1 + \cdots + X_n = \cfrac{1}{n} \displaystyle \sum_{i=1}^{n} X_n$$

を新しい確率変数とみなすと、$n$ を無限大に発散させたとき、$\bar{X}$ の確率分布は平均 $\mu$、分散 $\cfrac{\sigma^2}{n}$ の正規分布に収束する。

言い換えると、標本数 $n$ が十分に大きい場合、母集団の分布がどんな形であっても、標本平均 $\bar{X}$ の分布は正規分布へと近づいていく。

証明

平均・分散の証明

$\bar{X}$ の平均 $E(\bar{X}) = \mu$、分散 $V(\bar{X}) = \cfrac{\sigma^2}{n}$ であることは、一般的な期待値と分散の性質

$$E(aX) = aE(X)$$ $$E(X+Y) = E(X) + E(Y)$$ $$V(aX) = a^2 V(X)$$ $$V(X+Y) = V(X) + V(Y) + 2 \mathrm{Cov}(X, Y)$$

（ここで $\mathrm{Cov}$ は共分散であり、$X,Y$ が互いに独立ならゼロ）

を用いて、以下のように計算して示せる。

$$\begin{eqnarray} E(\bar{X}) &=& E \left( \displaystyle \cfrac{1}{n} \sum_{i=1}^{n} X_i \right) \\ &=& \cfrac{1}{n} \displaystyle \sum_{i=1}^{n} E(X_i) \\ &=& \cfrac{1}{n} \displaystyle \sum_{i=1}^{n} \mu \\ &=& \cfrac{1}{n} n\mu = \mu \end{eqnarray}$$ $$\begin{eqnarray} V(\bar{X}) &=& V \left( \displaystyle \cfrac{1}{n} \sum_{i=1}^{n} X_i \right) \\ &=& \cfrac{1}{n^2} \displaystyle \sum_{i=1}^{n} V(X_i) \\ &=& \cfrac{1}{n^2} \displaystyle \sum_{i=1}^{n} \sigma^2 \\ &=& \cfrac{1}{n^2} n \sigma^2 = \cfrac{\sigma^2}{n} \end{eqnarray}$$

正規分布に従うことの証明

$\bar{X}$ が正規分布に従う $\Longleftrightarrow$ $\bar{X}$ を標準化した確率変数 $Z$ が標準正規分布に従う

なので、$n \rightarrow \infty$ のときに $Z$ が標準正規分布に従うことを証明すれば良い。
それには、モーメント母関数（積率母関数）が標準正規分布のものに収束することを示せば十分。

前節で計算した平均値と分散の値を用いて $\bar{X}$ を標準化すると、

$$\begin{eqnarray} Z &\equiv& \cfrac{\bar{X}-\mu}{\sigma / \sqrt{n}} \\ &=& \cfrac{\sqrt{n}}{\sigma} \left( \cfrac{1}{n} \displaystyle \sum_{i=1}^{n} X_i - \mu \right) \\ &=& \cfrac{1}{\sigma \sqrt{n}} \left( \displaystyle \sum_{i=1}^{n} X_i - n \mu \right) \\ &=& \cfrac{1}{\sigma \sqrt{n}} \left( \displaystyle \sum_{i=1}^{n} (X_i - \mu) \right) \end{eqnarray}$$

ここで $X_i$ を標準化した確率変数

$$Z_i \equiv \cfrac{X_i-\mu}{\sigma}$$

を導入すれば、

$$Z = \cfrac{1}{\sqrt{n}} \displaystyle \sum_{i=1}^{n}Z_i$$

$Z$ のモーメント母関数は

$$\begin{eqnarray} M_Z(t) &=& E \left( \exp{(tZ)} \right) \\ &=& E \left( \exp{ \left( \cfrac{t}{\sqrt{n}} \displaystyle \sum_{i=1}^{n}Z_i \right) } \right) \\ &=& E \left( \displaystyle \prod_{i=1}^{n} \exp{ \left( \cfrac{t}{\sqrt{n}} Z_i \right) } \right) \\ &=& \displaystyle \prod_{i=1}^{n} E \left( \exp{ \left( \cfrac{t}{\sqrt{n}} Z_i \right) } \right) \\ &=& \displaystyle \prod_{i=1}^{n} M_{Z_i} \left (\cfrac{t}{\sqrt{n}} \right) \end{eqnarray}$$

よって、$Z$ のモーメント母関数は各変数 $Z_1, \cdots, Z_n$ のモーメント母関数の積で表現できる。

$Z_1, \cdots, Z_n$ は全て同じ確率分布に従うから、モーメント母関数も等しく、

$$M_Z(t) = \left( M_{Z_k} \left (\cfrac{t}{\sqrt{n}} \right) \right)^n$$

となる（$k$ は 1〜$n$ のどれでも良い）。

$M_{Z_k} \left (\cfrac{t}{\sqrt{n}} \right)$ をマクローリン展開すると、

$$M_{Z_k} \left (\cfrac{t}{\sqrt{n}} \right) = M_{X_k}(0) + \cfrac{M'_{Z_k}(0)}{1!}\left(\cfrac{t}{\sqrt{n}}\right) + \cfrac{M''_{Z_k}(0)}{2!}\left(\cfrac{t}{\sqrt{n}}\right)^2 + \cfrac{M'''_{Z_k}(0)}{3!}\left(\cfrac{t}{\sqrt{n}}\right)^3 + \cdots$$

ここで、

$$M_{Z_k} \left (\cfrac{t}{\sqrt{n}} \right) = \int_{-\infty}^{\infty} \exp{\left(\cfrac{tZ_k}{\sqrt{n}}\right)} f(Z_k) dZ_k$$

であるから、

$$M_{Z_k}(0) = \int_{-\infty}^{\infty} f(Z_k) dZ_k = 1$$

（任意の確率密度関数の全区間積分は1）

また、$Z_k$ は標準化されており、モーメント母関数の性質と分散・期待値の公式から、

$$M'_{Z_k}(0) = E(Z_k) = 0$$ $$M''_{Z_k}(0) = E(Z^2_k) = V(Z_k) + E(Z_k)^2 = 1+0^2=1$$

これらをマクローリン展開の式に代入し、$\cfrac{t}{\sqrt{n}}$ の3次以上の式を $O\left(\left(\cfrac{t}{\sqrt{n}}\right)^3\right)$ でまとめて表すと、

$$M_{Z_k} \left (\cfrac{t}{\sqrt{n}} \right) = 1 + \cfrac{t^2}{2n} + O\left(\left(\cfrac{t}{\sqrt{n}}\right)^3\right)$$

これを $M_Z(t)$ の式に代入すれば、

$$\begin{eqnarray} M_Z(t) &=& \left( M_{Z_k} \left( \cfrac{t}{\sqrt{n}} \right) \right)^n \\ &=& \left( 1 + \cfrac{t^2}{2n} + O\left(\left(\cfrac{t}{\sqrt{n}}\right)^3\right) \right)^n \\ &\longrightarrow& \exp{\left( \cfrac{t^2}{2} \right)} \qquad (n \rightarrow \infty) \end{eqnarray}$$

最後の式は標準正規分布のモーメント母関数に一致する（証明終了）。

実験

1. 色々な確率分布から $n$ 件の標本を繰り返し抽出して、標本平均の度数分布表を描画
2. 標本数 $n$ を大きくしていき、中心極限定理が成り立つ（= 標本数 $n$ が大きいほど標本平均の分布が正規分布に近づく）ことを確認

一様分布

ベータ分布

$\alpha = 9, \beta = 3$ の場合：

$\alpha = 2, \beta = 2$ の場合：

カイ二乗分布

自由度 $k = 2$ の場合：

自由度 $k = 5$ の場合：

一次関数に従う分布

ここでは、$0 \le x \le a$ の範囲で一次関数 $f(x) = cx$ 型の分布に従う母集団を考え、$a=1$ の場合を描画した。

確率密度関数 $f(x)$ の全区間積分（直角三角形の面積）は1になる必要があるので、三角形の高さ（$x=a$ における $f(x)$ の値）は $\cfrac{2}{a}$

したがって、この確率密度関数（直角三角形の斜辺）の傾き $c = \cfrac{2}{a^2}$ となるから、

$$f(x) = \begin{cases} \cfrac{2}{a^2} x & \qquad & (0 \le x \le a) \\ 0 & \qquad & (x \lt 0, a \lt x) \end{cases}$$

母集団の期待値 $\mu$ は

$$\begin{eqnarray} \mu &=& E(x) = \int_0^a x f(x) dx \\ &=& \cfrac{2}{a^2} \int_0^a x^2 dx = \cfrac{2}{a^2} \left[ \cfrac{x^3}{3} \right]_0^a \\ &=& \cfrac{2}{3} a \end{eqnarray}$$

同様に

$$\begin{eqnarray} E(x^2) &=& \int_0^a x^2 f(x) dx \\ &=& \cfrac{2}{a^2} \int_0^a x^3 dx = \cfrac{2}{a^2} \left[ \cfrac{x^4}{4} \right]_0^a \\ &=& \cfrac{1}{2} a^2 \end{eqnarray}$$

であるから、母集団の分散 $\sigma^2$ は

$$\begin{eqnarray} \sigma^2 &=& V(x) = E(x^2) - E(x)^2 \\ &=& \cfrac{1}{2} a^2 - \cfrac{4}{9} a^2 \\ &=& \cfrac{1}{18} a^2 \end{eqnarray}$$

二次関数に従う分布

ここでは、$-a \le x \le a$ の範囲で二次関数 $f(x) = cx^2$ 型の分布に従う母集団を考え、$a=1$ の場合を描画した。

確率密度関数 $f(x)$ の全区間積分は1になる必要があるので、

$$\begin{eqnarray} 1 &=& \int_{-a}^a f(x) dx = c \int_{-a}^a x^2 dx \\ &=& c \left[ \cfrac{x^3}{3} \right]_{-a}^a = \cfrac{2}{3} ca^3 \end{eqnarray}$$

よって分布関数の係数 $c = \cfrac{3}{2a^3}$ となるから、

$$f(x) = \begin{cases} \cfrac{3}{2a^3} x^2 & \qquad & (-a \le x \le a) \\ 0 & \qquad & (x \lt -a, a \lt x) \end{cases}$$

母集団の期待値 $\mu$ は

$$\begin{eqnarray} \mu &=& E(x) = \int_{-a}^a x f(x) dx \\ &=& \cfrac{3}{2a^3} \int_{-a}^a x^3 dx = \cfrac{3}{2a^3} \left[ \cfrac{x^4}{4} \right]_{-a}^a \\ &=& \cfrac{3}{8a^3} (a^4 - a^4) = 0 \end{eqnarray}$$

同様に

$$\begin{eqnarray} E(x^2) &=& \int_{-a}^a x^2 f(x) dx \\ &=& \cfrac{3}{2a^3} \int_{-a}^a x^4 dx = \cfrac{3}{2a^3} \left[ \cfrac{x^5}{5} \right]_{-a}^a \\ &=& \cfrac{3}{10a^3} (a^5 - (-a^5)) = \cfrac{3}{5} a^2 \end{eqnarray}$$

であるから、母集団の分散 $\sigma^2$ は

$$\begin{eqnarray} \sigma^2 &=& V(x) = E(x^2) - E(x)^2 \\ &=& \cfrac{3}{5} a^2 - 0^2 = \cfrac{3}{5} a^2 \end{eqnarray}$$

Appendix: 実験に使った Python コード

	import numpy as np
	import math
	from matplotlib import pyplot as plt
	from abc import ABC, abstractmethod
	def ln_factorial(n):
	"""階乗の log を計算するための補助関数"""
	ret = 0
	for i in range(1, n+1):
	ret += np.log(i)
	return ret
	def gauss(x, mu, sigma):
	"""x を引数として平均 mu, 標準偏差 sigma の正規分布の確率密度関数を計算"""
	ret = np.exp(- (x-mu)**2 * 0.5 / sigma**2)
	ret /= np.sqrt(2.0 * np.pi) * sigma
	return ret
	class Distribution(ABC):
	"""一般的な確率分布を表すクラス"""
	@abstractmethod
	def get_name(self):
	pass
	@abstractmethod
	def calc_dist_func_xy(self):
	"""確率密度関数描画のための x, y を返す"""
	pass
	@abstractmethod
	def get_average_of_samples(self, n):
	"""n 件の標本をランダムに抽出し、平均値を返す"""
	pass
	def get_mu(self):
	return self.__mu
	def get_sigma(self):
	return self.__sigma
	def set_mu(self, mu):
	self.__mu = mu
	def set_sigma(self, sigma):
	self.__sigma = sigma
	class UniformDist(Distribution):
	"""一様分布"""
	def __init__(self, a, b):
	self.set_mu((a + b) * 0.5)
	variance = (b-a)**2/12.0
	self.set_sigma(np.sqrt(variance))
	self.__a = a
	self.__b = b
	def get_name(self):
	return 'Uniform Distribution ($a={},b={}$)'.format(self.__a, self.__b)
	def calc_dist_func_xy(self):
	w = self.__b - self.__a
	x = np.linspace(self.__a-w, self.__b+w, 300)
	y = np.where((self.__a <= x) & (x <= self.__b), 1.0/w, 0)
	return x, y
	def get_average_of_samples(self, n):
	samples = np.random.rand(n) * (self.__b - self.__a) + self.__a
	return samples.mean()
	class BetaDist(Distribution):
	"""ベータ分布"""
	def __init__(self, a, b):
	self.set_mu(a/(a+b))
	variance = a*b/(a+b+1)/(a+b)**2
	self.set_sigma(np.sqrt(variance))
	self.__a = a
	self.__b = b
	def get_name(self):
	return r'Beta Distribution ($\alpha={}, \beta={}$)'.format(self.__a, self.__b)
	def calc_dist_func_xy(self):
	def beta(p):
	ret = np.zeros(p.shape)
	for i in range(len(p)):
	if 0 < p[i] < 1:
	log_beta = (self.__a-1) * np.log(p[i]) + (self.__b-1) * np.log(1-p[i]) + ln_factorial(self.__a+self.__b-1) - ln_factorial(self.__a-1) - ln_factorial(self.__b-1)
	ret[i] = np.exp(log_beta)
	return ret
	x = np.linspace(0, 1.0, 100)
	y = beta(x)
	return x, y
	def get_average_of_samples(self, n):
	samples = np.random.beta(self.__a, self.__b, n)
	return samples.mean()
	class ChiSquaredDist(Distribution):
	"""カイ二乗分布"""
	def __init__(self, k):
	self.set_mu(k)
	variance = 2 * k
	self.set_sigma(np.sqrt(variance))
	self.__k = k
	def get_name(self):
	return r'$\chi^2$ Distribution ($k={}$)'.format(self.__k)
	def calc_dist_func_xy(self):
	x = np.linspace(0, 20.0, 100)
	y = np.exp(-x*0.5) / (2**(self.__k*0.5)) / math.gamma(self.__k*0.5) * np.power(x, self.__k*0.5-1)
	return x, y
	def get_average_of_samples(self, n):
	samples = np.random.chisquare(self.__k, n)
	return samples.mean()
	class MyLinearDist(Distribution):
	"""一次関数型の分布"""
	def __init__(self, a):
	self.set_mu(a * 2.0 / 3.0)
	variance = a * a / 18.0
	self.set_sigma(np.sqrt(variance))
	self.__x_max = a
	self.__c = 2.0 / (a*a) # 直線の傾き
	self.__y_max = 2.0 / a # 分布内の y の最大値
	def get_name(self):
	return 'Self-Made Linear Distribution'
	def calc_dist_func_xy(self):
	x = np.linspace(0-self.__x_max*0.5, self.__x_max*1.5, 100)
	y = np.where((0 <= x) & (x <= self.__x_max), self.__c * x, 0)
	return x, y
	def get_average_of_samples(self, n):
	samples = []
	while len(samples) < n:
	x_tmp = np.random.rand() * self.__x_max
	y_tmp = np.random.rand() * self.__y_max
	if y_tmp <= self.__c * x_tmp:
	samples.append(x_tmp)
	return np.mean(samples)
	class MyQuadraticDist(Distribution):
	"""二次関数型の分布"""
	def __init__(self, a):
	self.set_mu(0)
	variance = a * a * 3.0 / 5.0
	self.set_sigma(np.sqrt(variance))
	self.__x_max = a
	self.__c = 3.0 / (2.0*a*a*a) # 二次関数の係数
	self.__y_max = 3.0 / (2.0*a) # 分布内の y の最大値
	def get_name(self):
	return 'Self-Made Quadratic Distribution'
	def calc_dist_func_xy(self):
	x = np.linspace(-self.__x_max*1.5, self.__x_max*1.5, 100)
	y = np.where((-self.__x_max <= x) & (x <= self.__x_max), x*x*self.__c, 0)
	return x, y
	def get_average_of_samples(self, n):
	samples = []
	while len(samples) < n:
	x_tmp = (np.random.rand()*2-1.0) * self.__x_max
	y_tmp = np.random.rand() * self.__y_max
	if y_tmp <= self.__c * x_tmp * x_tmp:
	samples.append(x_tmp)
	return np.mean(samples)
	def func(dist):
	"""
	与えられた確率分布について色々な標本サイズで標本平均を繰り返し計算し、
	標本平均の分布が正規分布に近づくことをグラフで確認する
	"""
	T = 100000 # n個の標本を抽出して平均を取る操作を何度繰り返すか
	n_samples = [2, 4, 8, 32, 128] # 標本サイズ
	plt.figure(figsize=(10, 6))
	plt.subplots_adjust(wspace=0.3, hspace=0.4)
	# 母集団の確率密度関数を描画
	plt.subplot(2, 3, 1)
	x, y = dist.calc_dist_func_xy()
	plt.title('Distribution Function')
	plt.grid()
	plt.fill_between(x, np.zeros(x.size), y, facecolor='red', alpha=0.3)
	plt.plot(x, y, color='black')
	# 標本を繰り返し抽出して平均値を記録・ヒストグラム描画
	for i in range(len(n_samples)):
	n = n_samples[i]
	sample_mean = np.empty(T, dtype='float')
	for t in range(T):
	sample_mean[t] = dist.get_average_of_samples(n)
	# ヒストグラム描画のための階級幅を計算
	mean_max = sample_mean.max()
	mean_min = sample_mean.min()
	bins = np.linspace(mean_min, mean_max, 50)
	# ヒストグラムを描画
	plt.subplot(2, 3, i+2)
	plt.hist(sample_mean, bins=bins, density=True)
	# 中心極限定理により収束が期待される正規分布を描画
	mu = dist.get_mu()
	sigma = dist.get_sigma() / np.sqrt(n)
	x_norm = np.linspace(mean_min-(mean_max-mean_min)*0.2, mean_max+(mean_max-mean_min)*0.2, 100)
	norm = gauss(x_norm, mu=mu, sigma=sigma)
	plt.title(r'$n_{{\rm sample}} = {}$'.format(n))
	plt.plot(x_norm, norm, lw=1.0, color='red', label=r'$N(\mu, \sigma^2/{})$'.format(n))
	plt.legend()
	plt.suptitle('Sample Mean of ' + dist.get_name())
	plt.show()
	func(UniformDist(0, 6.0))
	func(BetaDist(9, 3))
	func(BetaDist(2, 2))
	func(ChiSquaredDist(2))
	func(ChiSquaredDist(5))
	func(MyLinearDist(1.0))
	func(MyQuadraticDist(1.0))

view raw central-limit-theorem.py hosted with ❤ by GitHub