Up

正規分布、対数正規分布、中心極限定理とシミュレーション

多数の要因が働くと正規分布になるというように説明されることがある。これは、中心極限定理に基づくものである。正規分布は、要因の効果がそれらの和として現れる場合であるが、積として現れる場合は対数正規分布（Gelman et al., 2014, pp.576-577）になる。中心極限定理と正規分布、対数正規分布の関係を以下のシミュレーションで確認してみた。

中心極限定理は幾つかの形のものがあるが、例えば次のような説明がある（岡本、2009）。

確率変数が同じ分布に従う独立な確率変数であるとする。いま、

と置く。このとき、和はのとき標準正規分布に近づく。すなわち、

である。

上式の両辺の指数関数をとると

となる。

確率変数に、連続型として一様分布、離散型としてベルヌーイ分布をとったときの様子を以下のシミュレーションで調べた。

一様分分布の場合

n=100の場合のの分布を調べるシミュレーションプログラムを、リスト１のように作成した。

リスト１　一様分布の場合のシミュレーション

import scipy.stats as ss

import numpy as np

import matplotlib.pyplot as plt

n_simu = 1_000_000

n = 100

samples = []

for i in range(n_simu):

if i % 10000 == 0:

print(i, '/', n_simu)

x = ss.uniform.rvs(size=n)

x_n = (x - 1/2) * 2 * (3**0.5) / (n**0.5)

sum_x_n = np.sum(x_n)

samples.append(sum_x_n)

min_x = np.min(samples)

max_x = np.max(samples)

Bins = np.linspace(min_x, max_x, 100)

xcoord = Bins + (max_x - min_x) / (100*2)

ycoord = ss.Normal(mu=0, sigma=1).pdf(xcoord)

plt.hist(samples, bins = Bins, density=True)

plt.plot(xcoord, ycoord)

plt.show()

次のコードにより、n個の一様乱数を取り、に変換している。

x = ss.uniform.rvs(size=n)

x_n = (x - 1/2) * 2 * (3**0.5) / (n**0.5)

次のコードにより、変換したものの和を取り、リストsamplesに格納している。

sum_x_n = np.sum(x_n)

samples.append(sum_x_n)

次のコードにより、samplesのヒストグラムと、比較のための標準正規分布のグラフを表示している。

Bins = np.linspace(min_x, max_x, 100)

xcoord = Bins + (max_x - min_x) / (100*2)

ycoord = ss.Normal(mu=0, sigma=1).pdf(xcoord)

plt.hist(samples, bins = Bins, density=True)

plt.plot(xcoord, ycoord)

plt.show()

グラフは、図１のように表示される。

図１

N=100の場合のグラフである。ヒストグラムは、標準正規分布のグラフによく一致している。

次に、が対数正規分布に近づくことを確認する。スクリプトはリスト２のように作成した。

リスト２　一様分布を用いて積が対数正規分布に近づくことを見るスクリプト

import scipy.stats as ss

import numpy as np

import matplotlib.pyplot as plt

n_simu = 1_000_000

n = 100

samples = []

for i in range(n_simu):

if i % 10_000 == 0:

print(i, '/', n_simu)

x =ss.uniform.rvs(size=n)

x_n = (x - 1/2) * 2 * (3**0.5) / (n**0.5)

y_n = np.exp(x_n)

prod_y_n = np.prod(y_n)

samples.append(prod_y_n)

min_x = np.min(samples)

max_x = np.max(samples)

Bins = np.linspace(min_x, max_x, 100)

xcoord = Bins + (max_x - min_x) / (100*2)

ycoord = ss.lognorm.pdf(xcoord, 1)

plt.hist(samples, bins = Bins, density=True)

plt.plot(xcoord, ycoord)

plt.show()

samples = np.array(samples)

sub_samples = samples[np.argsort(samples)]

sub_samples = sub_samples[:-int(n_simu/100)]

sub_max_x = np.max(sub_samples)

sub_Bins = np.linspace(min_x, sub_max_x, 100)

sub_xcoord = sub_Bins + (sub_max_x - min_x) / (100*2)

sub_ycoord = ss.lognorm.pdf(sub_xcoord, 1) / 0.99

plt.hist(samples, bins = sub_Bins, density=True)

plt.plot(sub_xcoord, sub_ycoord)

plt.show()

を一様分布としたときの積の計算を次のコードで行っている。

x =ss.uniform.rvs(size=n)

x_n = (x - 1/2) * 2 * (3**0.5) / (n**0.5)

y_n = np.exp(x_n)

prod_y_n = np.prod(y_n)

samples.append(prod_y_n)

上のコードによるシミュレーションの結果を、次のコードで表示すると図２のような結果を得る。

xcoord = Bins + (max_x - min_x) / (100*2)

ycoord = ss.lognorm.pdf(xcoord, 1)

plt.hist(samples, bins = Bins, density=True)

plt.plot(xcoord, ycoord)

plt.show()

図２

右裾が長く、分布の大半が左端によっているので、サンプルの１%を右裾から削って描画してみる。

図２のグラフが得られる。

図３

ヒストグラムは、橙色の曲線で描かれている対数正規分布とよく合っていることがわかる。

ベルヌーイ分布の場合

をベルヌーイ分布Bernoulli(0.5)に従うとしたときのシミュレーションのスクリプトをリスト３のように用意した。

リスト３　ベルヌーイ分布の場合の和のシミュレーション

import scipy.stats as ss

import numpy as np

import matplotlib.pyplot as plt

n_simu = 1_000_000

n = 1000

samples = []

for i in range(n_simu):

if i % 10_000 == 0:

print(i, '/', n_simu)

x = ss.bernoulli.rvs(0.5, size=n)

x_n = (x - 1/2) * 2 / (n**0.5)

sum_x_n = np.sum(x_n)

samples.append(sum_x_n)

min_x = np.min(samples)

max_x = np.max(samples)

Bins = np.linspace(min_x, max_x, 20)

xcoord = Bins + (max_x - min_x) / (20*2)

ycoord = ss.Normal(mu=0, sigma=1).pdf(xcoord)

plt.hist(samples, bins = Bins, density=True)

plt.plot(xcoord, ycoord)

plt.show()

和を次のコードで作成してリストsamplesに格納している。

x = ss.bernoulli.rvs(0.5, size=n)

x_n = (x - 1/2) * 2 / (n**0.5)

sum_x_n = np.sum(x_n)

samples.append(sum_x_n)

和の個数nはn=1000である。

リストsamplesのヒストグラムを、次のコードで描くと図４のようになる。

plt.hist(samples, bins = Bins, density=True)

図４

橙色の曲線は、標準正規分布である。ヒストグラムが標準正規分布をよく近似していることがわかる。

次に、積のシミュレーションを行う。

リスト４のスクリプトを作成した。

リスト４　積のシミュレーション

import scipy.stats as ss

import numpy as np

import matplotlib.pyplot as plt

n_simu = 1_000_000

n = 1000

samples = []

for i in range(n_simu):

if i % 10_000 == 0:

print(i, '/', n_simu)

x = ss.bernoulli.rvs(0.5, size=n)

x_n = (x - 1/2) * 2 / (n**0.5)

y_n = np.exp(x_n)

prod_y_n = np.prod(y_n)

samples.append(prod_y_n)

min_x = np.min(samples)

max_x = np.max(samples)

Bins = np.linspace(min_x, max_x, 20)

xcoord = Bins + (max_x - min_x) / (20*2)

ycoord = ss.lognorm.pdf(xcoord, 1)

plt.hist(samples, bins = Bins, density=True)

plt.plot(xcoord, ycoord)

plt.show()

samples = np.array(samples)

sub_samples = samples[np.argsort(samples)]

sub_samples = sub_samples[:-int(n_simu/100)]

sub_max_x = np.max(sub_samples)

sub_Bins = np.linspace(min_x, sub_max_x, 20)

sub_xcoord = sub_Bins + (sub_max_x - min_x) / (20*2)

sub_ycoord = ss.lognorm.pdf(sub_xcoord, 1) / 0.99

ck_Bins = np.linspace(min_x, sub_max_x, 1000)

ck_ycoord = ss.lognorm.pdf(ck_Bins, 1) / 0.99

plt.hist(samples, bins = sub_Bins, density=True, alpha=0.5)

plt.plot(sub_xcoord, sub_ycoord, label='n_grid=20')

plt.plot(ck_Bins, ck_ycoord, label='n_grid=1000')

plt.legend()

plt.show()

積の計算と、計算結果のリストsamplesへの格納を次のコードで行っている。

x = ss.bernoulli.rvs(0.5, size=n)

x_n = (x - 1/2) * 2 / (n**0.5)

y_n = np.exp(x_n)

prod_y_n = np.prod(y_n)

samples.append(prod_y_n)

作成したリストsamplesのヒストグラム描画のコードを以下のように用意した。

xcoord = Bins + (max_x - min_x) / (20*2)

ycoord = ss.lognorm.pdf(xcoord, 1)

plt.hist(samples, bins = Bins, density=True)

plt.plot(xcoord, ycoord)

plt.show()

上のコードにより、図５のようなグラフが描かれる。

図５

右に裾が長く伸びているので、右側１%を除いて描画するコードを以下のように用意した。

samples = np.array(samples)

sub_samples = samples[np.argsort(samples)]

sub_samples = sub_samples[:-int(n_simu/100)]

sub_max_x = np.max(sub_samples)

sub_Bins = np.linspace(min_x, sub_max_x, 20)

sub_xcoord = sub_Bins + (sub_max_x - min_x) / (20*2)

sub_ycoord = ss.lognorm.pdf(sub_xcoord, 1) / 0.99

ck_Bins = np.linspace(min_x, sub_max_x, 1000)

ck_ycoord = ss.lognorm.pdf(ck_Bins, 1) / 0.99

plt.hist(samples, bins = sub_Bins, density=True, alpha=0.5)

plt.plot(sub_xcoord, sub_ycoord, label='n_grid=20')

plt.plot(ck_Bins, ck_ycoord, label='n_grid=1000')

plt.legend()

plt.show()

上のコードにより、図６のようなグラフが描かれる。

図６

橙色の折れ線は、ヒストグラムのビンの中点における対数正規分布の確率密度関数の値を結んだものである。ヒストグラムでは、左端の０に近いところでの頻度の低いところの影響により、ヒストグラムと折れ線との乖離が顕著である。これをグリッドを細かくして計算した緑色の曲線では、この左のビン内でのヒストグラムと対数正規分分布の曲線との関係を見ることができる。

参考文献

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D> B., Vehtari, A, & Rubin, D. B. (2014). Bayesian data analysis, 3rd ed. CRC Press (free download)

岡本安晴（2009）　統計学を学ぶための数学入門［下］．培風館

Home