Pythonで〇×ゲームのAIを一から作成する その231 大数の法則と中心極限定理の関係

目次と前回の記事

Python のバージョンとこれまでに作成したモジュール

本記事のプログラムは Python のバージョン 3.13 で実行しています。また、numpy のバージョンは 2.3.5 です。

リンク	説明
marubatsu.py	Marubatsu、Marubatsu_GUI クラスの定義
ai.py	AI に関する関数
mbtest.py	テストに関する関数
util.py	ユーティリティ関数の定義
tree.py	ゲーム木に関する Node、Mbtree クラスなどの定義
gui.py	GUI に関する処理を行う基底クラスとなる GUI クラスの定義

AI の一覧とこれまでに作成したデータファイルについては、下記の記事を参照して下さい。

今回の記事の内容

次回の記事で説明する予定の原始モンテカルロ法を用いた 〇× ゲームの AI の作成には必要ありませんが、大数の法則についてこれまで詳しく説明したので、大数の法則に密接に関連する、統計学の中で最も重要な定理の一つである中心極限定理について説明します。

中心極限定理

大数の法則は母集団から無作為復元抽出された標本の標本平均が、標本サイズを大きくすると母平均に近づく（確率収束する）というものですが、標本サイズを大きくした際の標本平均の分布がどのように変化するかについては言及されていません。そのことを具体的に表す定理が中心極限定理です。

中心極限定理の定義を要約すると以下のようになります。

有限な値の母平均 $μ$ と母分散 $σ^2$ を持つ母集団の確率分布を表す確率変数を $X$ とし、その母集団から $n$ 個の要素を無作為復元抽出を行うものとする。作成された標本サイズが $n$ の標本の標本平均を表す確率変数を $\bar{X_n}$ とすると、$n$ を十分大きくすると $\bar{X_n}$ が従う確率分布は期待値が $μ$、分散が $\frac{σ ^ 2}{n}$ の正規分布で近似できる。

$\bar{X_n}$ の期待値 $E[\bar{X_n}]$ が $μ$ であることは以前の記事で、分散 $V[\bar{X_n}]$ が $\frac{σ ^ 2}{n}$ であることは以前の記事で証明済です。中心極限定理ではそのことに加えて、標本平均の 確率分布の形状 が正規分布に近づいていくことを示し、このような性質の事を分布収束と呼びます。

上記で「$n$ を十分大きくすると」という曖昧な表現になっているのは、元の確率分布によってどれくらい $n$ 大きくなれば正規分布に近似できるかが異なるからです。なお、一般的には $n$ が 30 以上になると正規分布で近似しても良いと言われることが多いようです。

間違えやすいですが、中心極限定理によって正規分布に近づいていくのは元の $X$ の確率分布ではなく、標本平均 $\bar{X_n}$ の確率分布であるという点に注意して下さい。

当然ですが標本サイズをいくら大きくしたからといって、元の確率分布が変化することはありません。例えばサイコロを何回振ってその出目を記録（無作為復元抽出）しても、サイコロを振った時の出目の確率は変化しません。

中心極限定理が非常に便利な点は、有限な期待値と分散を持つという条件を満たせば元の確率分布がどのような分布であっても成り立つという点にあります。どのように便利であるかの具体例は今回の記事の後半で示します。

なお、中心極限定義の証明は非常に難しいので、今回の記事では中心極限定理が様々な確率分布に対して成り立つことを確認することにします。

大数の法則は確率分布が有限な分散を持つことを条件としていませんが、中心極限定理では有限な分散を持つことを条件とする点が異なります。

ただし、有限な分散を持たない場合でも大数の法則が満たされることから想像できるように、確率分布が有限の期待値を持つという条件を満たせば正規分布とは異なる特定の確率分布に分布収束します。また、有限な分散を持つという条件を必要としない場合の中心極限定理の事を一般化中心極限定理と呼びます。

参考までに Wikipedia の 中心極限定理の項目を下記に示します。

標本平均を表す確率変数の分散の性質

標本サイズ $n$ を大きくすると標本平均を表す確率変数 $\bar{X_n}$ のばらつきを表す分散がどのように変化するかを説明します。以前の記事で証明したように、標本サイズが $n$ の標本平均の分散 $V[\bar{X_n}]$ は下記の式のように母分散 $σ^2$ を $n$ で割った $\frac{σ^2}{n}$ になります。

$V[\bar{X_n}] = \frac{V[X]}{n} = \frac{σ^2}{n}$

中心極限定理の前提条件から $σ^2$ は有限な値なので、標本サイズ $n$ を大きくしていくと $V[\bar{X_n}]$ の極限値は下記の式のように 0 になります。

$\lim_{n \to ∞} V[\bar{X_n}] = \lim_{n \to ∞} \frac{σ^2}{n} = 0$

以前の記事で説明したように、分散が 0 であるということはばらつきが全くない（すべての実現値が平均と同じになる）ことを表します。標本平均 $E[\bar{X_n}]$ の期待値は常に母平均である $μ$ に一致するので、この式は大数の法則のメカニズムを下記のように説明してくれます。

• $n$ を大きくすると標本平均のばらつきを表す $V[\bar{X_n}]$ がどんどん 0 に向かって小さくなる
• その結果、標本平均 $\bar{X_n}$ の実現値が母平均 $μ$ のごく近くに集中するようになる
• 従って、標本サイズ $n$ を十分大きく取れば 1 回の標本抽出で得られる標本の標本平均はほぼ確実に母平均 $μ$ の付近の値を取る（大数の法則）

標本平均の分散を計算する関数の定義

上記が成り立つことを、プログラムによるシミュレーションで確認してみることにします。そのために、下記の計算を行う sampling_and_analyze という関数を定義することにします。標本抽出（sampling）と分析（analyze）を行うのでそのような名称にしました。

• 指定した確率分布の母分散 $σ^2$ を計算する
• その確率分布から無作為復元抽出によって標本サイズが size の標本を num 個作成する
• 作成した num 個の標本の標本平均を計算する
• num 個の標本平均の分散を計算する
• 母分散と $\frac{σ^2}{n}$ を表す「標本平均の分散 ÷ size」 を比較する

下記は、sampling_and_analyze の仮引数の一覧です。X と P は 以前の記事で定義した create_pd が作成する確率分布のデータと同じものです。なお、用語が似ているので間違えやすいですが、「標本サイズ」は 1 つの標本の中の要素の数を、「標本数」は作成した標本の数を表す点に注意して下さい。

仮引数	意味
X	確率変数 $X$ がとりうる値の一覧を表す 1 次元の ndarray
P	X の各要素の値が発生する確率
size	作成する標本の標本サイズ
num	作成する標本数

母平均と母分散の計算方法

最初に X と P から母平均 $μ = E[X]$ と母分散 $σ^2 = V[X]$ を計算する必要があるので、その計算方法について説明します。

母平均は numpy の np.average という関数で計算することができます。np.average は指定した ndarray の要素の平均（average）を計算する関数ですが、その際に仮引数 weights に各要素の重み（weight）を代入することで、「各要素に重みを乗算した値の合計」÷「重みの合計」という、重みを考慮した平均値を計算します。

np.average の詳細は下記のリンク先を参照して下さい。

下記のプログラムの 1 行目は 1 と 2 の平均である 1.5 を計算しますが、2 行目は 1 × 1 と 2 × 3 の合計である 7 を、重みの合計である 4 で割った 1.75 を計算します。

import numpy as np

print(np.average([1, 2]))
print(np.average([1, 2], weights=[1, 3]))

実行結果

1.5
1.75

確率分布の期待値は、確率分布の各要素の値とその確率を乗算した値の合計で計算でき、確率を表す P の合計は 1 になるため、np.average の仮引数 weights に P を代入することで X と P で表される確率分布の期待値を計算することができます。下記は偏りのないコインを投げた時の表を 0、裏を 1 とした場合の確率分布の期待値を計算するプログラムです。実行結果から 0.5 という正しい計算が行われることが確認できます。

X = np.array([0, 1])
P = np.array([0.5, 0.5])
print(np.average(X, weights=P))

以前の記事で説明したように、分散は下記の公式で計算することができます。このうちの $E[X]$ は上記のプログラムで計算できます。

$V[X] = E[X^2] - E[X]^2$

$V[X] = E[(X- E[X])^2]$ で計算することもできます。興味がある方はこの式で計算してみて下さい。

$E[X^2]$ は先ほどの X の代わりに X の各要素の 2 乗を計算する X ** 2 を記述した下記のプログラムで計算することができます。X * X は [0 * 0, 1 * 1] = [0, 1] なのでその期待値は上記と同じ 0.5 になり、実行結果から正しい計算が行われることが確認できます。

print(np.average(X ** 2, weights=P))

実行結果

0.5

従って、X と P で表される確率分布の分散は下記のプログラムで計算することができます。実行結果から $0.5 - 0.5^2 = 0.25$ が正しく計算されることが確認できます。

print(np.average(X ** 2, weights=P) - np.average(X, weights=P) ** 2)

実行結果

0.25

標本平均をまとめて計算する方法

次に、無作為復元抽出によって標本サイズが size の標本を num 個作成し、その標本平均の分散を計算する方法を説明します。

前回の記事で説明したように、無作為復元抽出による標本サイズが size の標本は下記のプログラムで計算することができます。

np.random.choice(X, p=P, size=size) 

上記のプログラムを num 回繰り返すことでこの標本を num 個作成することもできますが、下記のプログラムのように仮引数 size に (size, num) という tuple を代入することで、num 個の標本を表す (size, num) という shape を持つ 2 次元の ndarray を 1 回の np.random.choice でまとめて計算することもできます。

size = 5
num = 10
sample = np.random.choice(X, p=P, size=(size, num))
print(sample)

実行結果

[[0 1 0 0 0 0 0 0 0 1]
 [0 1 0 1 0 1 0 1 0 0]
 [1 1 0 1 1 1 0 0 0 1]
 [0 0 1 1 1 0 1 1 0 1]
 [1 1 0 1 0 0 1 0 0 0]]

以前の記事で説明したように、多次元の ndarray のそれぞれのインデックスを axis と呼び、先頭から axis 0、axis 1、・・・ のように番号で区別して表記します。上記の実行結果では縦方向が axis 0 に、横方向が axis 1 に対応します。

上記で計算される 2 次元の ndarray は最初のインデックスを表す axis 0 が size 個ある標本の要素のうちの 何番の要素 であるかを表し、次のインデックスを表す axis 1 が num 個ある標本のうちの 何番の標本 であるかを表します。例えば sample[0][0] は 0 番の標本の 0 番の要素を、sample[1][5] は 5 番の標本の 1 番の要素であることを表します。なお、要素の番号は 0 から数えるので最初の要素は 0 番である点に注意して下さい。従って、上記の 2 次元の ndarray のデータから、標本の各要素を表す axis 0 を対象に平均を計算することで、10 個の標本の平均をまとめて計算することができます。

以前の記事で説明した np.flip と同様に、numpy の多くの関数は仮引数 axis に値を代入することで、特定の axis を対象とした計算を行うことができます。np.average の場合は下記のプログラムのように仮引数 axis に 0 を代入することで axis 0 を対象とした平均を計算し、実行結果には 10 個の標本の標本平均の一覧が計算された ndarray が表示されます。先程の実行結果の各列（縦方向に並んだ 5 つの要素）の平均がそれぞれ計算されていることを確認して下さい。なお、要素の平均（mean）を計算すれば良いので、先程と異なり重みを表す仮引数 weights に値を代入する必要はありません。

mean = np.average(sample, axis=0)
print(mean)

実行結果

[0.4 0.8 0.2 0.8 0.4 0.4 0.4 0.4 0.  0.6]

np.average と似た関数に np.mean があります。np.mean は重みを考慮した計算を行うことができない点や、マスク付き ndarray に対応している点などが異なります。

np.mean の詳細は下記のリンク先を参照して下さい。

複数の標本をまとめて計算したほうが処理速度が高速になります。下記は size と num を 100 とした場合にまとめて計算を行う場合と、100 回の繰り返し処理によって 100 個の標本を作成する場合の処理時間を計測するプログラムです。実行結果から、まとめて計算を行うほうが処理速度が 10 倍以上になることが確認できました。

size = 100
num = 100

%timeit np.random.choice(X, p=P, size=(size, num))
%timeit [np.random.choice(X, p=P, size=size) for _ in range(num)]

実行結果

176 μs ± 461 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
2.26 ms ± 32.3 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

分散の計算方法

1 次元の ndarray の分散は下記のプログラムのように np.var という関数で計算できます。

print(np.var(mean))

実行結果

0.05600000000000001

np.var の詳細は下記のリンク先を参照して下さい。

以前の記事で説明した不偏分散は np.var の仮引数 ddof に 1 を代入することで計算できます。

sampling_and_analyze の定義

下記は、上記で説明した方法で「母平均と母分散の計算」、「複数個の標本の無作為復元抽出」、「それぞれの標本の標本平均とその分散の計算」を行う sampling_and_analyze を定義するプログラムです。

• 1 行目：先程説明した仮引数を持つ関数を定義する
• 2 行目：標本サイズを表示する
• 3 ~ 6 行目：母平均と母分散を計算して表示する
• 7、8 行目：標本サイズが size の標本を num 個抽出し、num 個の標本平均を計算する
• 9 行目：標本平均の分散を計算する
• 10 ~ 12 行目：標本平均の分散、母分散 ÷ 標本サイズ、その 2 つの誤差を表示する。なお、誤差は絶対値を計算する np.abs を利用して正の値として表示するようにした

 1  def sampling_and_analyze(X, P, size, num):
 2      print(f"標本サイズ = {size}")
 3      E = np.average(X, weights=P)
 4      V = np.average(X ** 2, weights=P) - E ** 2
 5      print(f"母平均           = {E:.7f}")
 6      print(f"母分散           = {V:.7f}")
 7      sample = np.random.choice(X, p=P, size=(size, num)) 
 8      mean = np.mean(sample, axis=0)  
 9      var = np.var(mean)
10      print(f"標本平均の分散    = {var:.7f}") 
11      print(f"母分散÷標本サイズ = {V/size:.7f}") 
12      print(f"誤差             = {np.abs(var - V/size):.7f}") 

行番号のないプログラム

def sampling_and_analyze(X, P, size, num):
    print(f"標本サイズ = {size}")
    E = np.average(X, weights=P)
    V = np.average(X ** 2, weights=P) - E ** 2
    print(f"母平均           = {E:.7f}")
    print(f"母分散           = {V:.7f}")

    sample = np.random.choice(X, p=P, size=(size, num)) 
    mean = np.mean(sample, axis=0)  
    var = np.var(mean)
    print(f"標本平均の分散    = {var:.7f}") 
    print(f"母分散÷標本サイズ = {V/size:.7f}") 
    print(f"誤差             = {np.abs(var - V/size):.7f}") 

np.abs の詳細は下記のリンク先を参照して下さい。

分割処理を行うようにする改良

上記のプログラムは 7 行目で shape が (size, num) の 2 次元の ndarray を作成します。その ndarray の要素の数は size * num なので、例えば標本サイズを 1 万、標本数を 1 万 とした場合は 1 億もの要素を持つ ndarray が生成されてしまい、前回の記事で説明したのと同様にコンピューターのメモリ不足が発生する可能性があります。

そこで前回の記事と同様に、一度に作成する ndarray の要素の数を制限して分割処理を行うように改良を行うことにします。

抽出する標本の要素の数に対する分割処理

前回の記事では一度の処理で 1 つの標本だけを抽出していたので、一度で抽出する標本の要素の数を 100 万に制限しました。今回の記事では複数の標本をまとめて抽出する点が異なりますが、最初は標本サイズだけに注目し、前回の記事と同様に一度に抽出する標本の要素の数を特定の数で制限するという分割処理を行うことにします。

下記は一度の処理で抽出する標本の要素（element）の数の上限を表す変数 maxe の値に応じた分割処理を行うプログラムです。また、その際に標本サイズ size = 8、標本数 num = 3、maxe = 5 を具体例として説明します。確率分布は先ほどの偏りのないコインを投げた場合のものを利用しました。

sizeleft や s を計算する処理は前回の記事のプログラムと同様ですが、この修正によって標本の要素を分割して抽出するようになるため、標本平均を np.average で計算することができなくなります。そのため、分割した標本の要素の合計を計算し、最後に標本サイズで割り算することで標本平均を計算する必要があります。

• 4 行目：抽出する必要がある残りの標本の要素数を表す変数 sizeleft を size で初期化する
• 5 行目：num 個の各標本の要素の合計を計算する変数を以前の記事で説明した np.zeros を利用して、要素数が num ですべての値が 0 である 1 次元の ndarray で初期化する
• 6 ~ 15 行目：sizeleft が 0 より大きい場合に繰り返し処理を行うことで、必要な標本サイズの標本を分割処理で抽出する
• 7、8 行目：一度に抽出する標本の標本サイズを表す s を計算し、その標本サイズで num 個の標本をまとめて抽出する
• 9 ~ 13 行目：計算の途中経過を確認できるように、抽出したデータと抽出したデータ内のそれぞれの標本の合計を表示する。次の繰り返し処理が区別できるように 13 行目で空行を表示するようにした。なお、抽出したデータ内の標本の合計は、先程標本平均を計算した np.average の代わりに np.sum を記述することで計算できる
• 14 行目：total に抽出したそれぞれの標本の合計を加算する
• 15 行目：s 個の要素を抽出したので、sizeleft からその値を引き算する
• 16 行目：繰り返し処理が終了した時点で total に各標本の要素の合計が計算されるのでそれを表示する
• 17、18 行目：各標本の要素の合計を標本サイズで割り算することで各標本の平均を計算して表示する
• 19、20 行目：np.var で標本平均の分散を計算して表示する

 1  size = 8
 2  num = 3
 3  maxe = 5
 4  sizeleft = size
 5  total = np.zeros(num)
 6  while sizeleft > 0:
 7      s = min(sizeleft, maxe)
 8      sample = np.random.choice(X, p=P, size=(s, num)) 
 9      print("sample")
10      print(f"{sample}")
11      print("sum")
12      print(f"{np.sum(sample, axis=0)}")
13      print()
14      total += np.sum(sample, axis=0)  
15      sizeleft -= s
16  print(f"total = {total}")
17  mean = total / size
18  print(f"mean  = {mean}")
19  var = np.var(mean)
20  print(f"var   = {var}")

行番号のないプログラム

size = 8
num = 3
maxe = 5
sizeleft = size
total = np.zeros(num)
while sizeleft > 0:
    s = min(sizeleft, maxe)
    sample = np.random.choice(X, p=P, size=(s, num)) 
    print("sample")
    print(f"{sample}")
    print("sum")
    print(f"{np.sum(sample, axis=0)}")
    print()
    total += np.sum(sample, axis=0)  
    sizeleft -= s
print(f"total = {total}")
mean = total / size
print(f"mean  = {mean}")
var = np.var(mean)
print(f"var   = {var}")

実行結果

sample
[[1 1 0]
 [1 0 1]
 [0 1 0]
 [1 1 1]
 [0 0 0]]
sum
[3 3 2]

sample
[[1 1 1]
 [1 1 0]
 [0 1 0]]
sum
[2 3 1]

total = [5. 6. 3.]
mean  = [0.625 0.75  0.375]
var   = 0.024305555555555552

s = min(sizeleft, maxe) = min(8, 5) = 5 なので、1 回目の繰り返しの 8 行目の処理では「5 つの要素を持つ 3 つの標本」を表す (5, 3) の shape の ndarray が計算され、その値が上記の実行結果の最初の sample として表示されます。また、3 つの標本の要素の合計がその下の sum に正しく計算されていることが確認できます。

上記の抽出を行った結果、抽出する必要がある残りの標本の要素数は sizeleft = 8 - 5 = 3 になります。s = min(sizeleft, maxe) = min(3, 5) = 3 なので、8 行目の処理では残りの「3 つの要素を持つ 3 つの標本」を表す (3, 3) の shape の ndarray が計算され、その値が上記の実行結果の 2 つ目の sample として表示されます。また、3 つの標本の要素の合計がその下の sum に正しく計算されていることが確認できます。

繰り返し処理の終了後の total には実行結果に表示される 2 つの sum の要素ごとの合計が計算されていることが確認でき、その下の mean には total の各要素を標本サイズを表す size = 8 で割った標本平均が計算されていることも確認できます。

実行結果の最後には標本平均の一覧から np.var を利用して標本平均の分散が計算されて表示されることが確認できます。

抽出する標本に関する分割処理

次に、一度に抽出する標本の数を制限することで、一度に np.random.choice によって計算される ndarray の要素数の上限が maxe になるように分割処理を行うことにします。

具体的には、一度に抽出する標本の数の最大値を表す maxn を下記の式で計算します。

maxs = min(size, maxe)
maxn = min(num, maxe // maxs)

最初の maxs は一度で抽出する標本の要素の数の最大値を表します。先ほどの処理では s = min(sizeleft, maxe) という式で計算しましたが、sizeleft の初期値は size で、繰り返し処理を行うたびに sizeleft の値は減っていくので、min(size, maxe) が s の最大値になります。

np.random.choice が計算する 2 次元の ndarray の要素の数は ndarray の shape が (s, n) の場合は s * n になります。そのため s * n の上限を maxe とするためには、s の最大値が maxs の場合は n の最大値を maxe // maxs 以下とする必要があります。なお、割り算の商を表す // で計算を行う理由は、n が整数である必要があるからです。

上記から maxn = min(num, maxe // maxs) という計算式で一度に抽出する標本の数の上限を計算することで、np.random.choice で計算する 2 次元の ndarray の要素の数の上限が maxe になります。

わかりづらいと思いますので具体例を挙げて説明します。標本サイズ size = 10000、標本数 num = 500、要素の最大値 maxe = 1000000（百万）とした場合の maxs と maxn は下記のように計算されます。

• maxs = min(size, maxe) = min(10000, 1000000) = 10000 となる
• maxn = min(num, maxe // maxs) = min(500, 1000000 // 10000) = min(500, 100) = 100 となる

shape を (size, num) とした 2 次元の ndarray の要素数は 10000 * 500 = 500 万 ですが、上記から (maxs, maxn) = (10000, 100) の shape を持つ ndarray の要素数は上限を表す maxe と同じ 100 万になります。

size = 100、num = 50、maxe = 1000000（百万）の場合の maxs と maxn は下記のようになります。

• maxs = min(100, 1000000) = 100 となる
• maxn = min(50, 1000000 // 100) = min(50, 10000) = 50 となる

この場合は (size, num) とした 2 次元の ndarray の要素数は 100 * 50 = 5000 で 100 万以下なので、(maxs, maxn) は (size, num) と完全に一致します。

下記はそのような処理を行うプログラムです。

• 5 行目：抽出する必要がある残りの標本数を表す変数 numleft を num で初期化する
• 6、7 行目：先程の計算式に基づき maxs と maxn を計算する
• 8 行目：各標本の標本平均を記録する変数を空の list で初期化する
• 9 ~ 27 行目：numleft が 0 より大きい間繰り返し処理を行うことで、必要な数の標本が抽出されるまで繰り返し処理を行う
• 10 ~ 22 行目：先程の処理と同様だが、11 行目で numleft と maxn の最小値を計算して n に代入し、15 行目の num を n とする点が異なる
• 24 行目：各標本の標本平均を表す list の要素に、上記で計算した標本平均の一覧を表すデータを + 演算子で結合して追加する。下記のノートで説明するが、その際に ndarray の tolist メソッドを利用して 1 次元の ndarray を list に変換する必要がある
• 25 行目：確認のためにこれまでに計算した標本平均の一覧を表示する
• 26 行目：n 個の標本を抽出したので numleft から n を引き算する
• 27 行目：次の標本を抽出することを明確にするために print で空行を表示する
• 28、29 行目：np.var で標本平均の分散を計算する。なお、np.var に限らず多くの numpy の関数は list を実引数に記述しても正しい計算を行うことができる

 1  size = 8
 2  num = 3
 3  maxe = 5
 4  
 5  numleft = num
 6  maxs = min(size, maxe)
 7  maxn = min(num, maxe // maxs)
 8  mean = []
 9  while numleft > 0:
10      sizeleft = size
11      n = min(numleft, maxn)
12      total = np.zeros(n)
13      while sizeleft > 0:
14          s = min(sizeleft, maxe)
15          sample = np.random.choice(X, p=P, size=(s, n)) 
16          print("sample")
17          print(f"{sample}")
18          print("sum")
19          print(f"{np.sum(sample, axis=0)}")
20          print()
21          total += np.sum(sample, axis=0)  
22          sizeleft -= s
23      print(f"total  = {total}")
24      mean += (total / size).tolist()
25      print(f"mean   = {mean}")
26      numleft -= n
27      print()
28  var = np.var(mean)
29  print(f"var    = {var}")

行番号のないプログラム

size = 8
num = 3
maxe = 5

numleft = num
maxs = min(size, maxe)
maxn = min(num, maxe // maxs)
mean = []
while numleft > 0:
    sizeleft = size
    n = min(numleft, maxn)
    total = np.zeros(n)
    while sizeleft > 0:
        s = min(sizeleft, maxe)
        sample = np.random.choice(X, p=P, size=(s, n)) 
        print("sample")
        print(f"{sample}")
        print("sum")
        print(f"{np.sum(sample, axis=0)}")
        print()
        total += np.sum(sample, axis=0)  
        sizeleft -= s
    print(f"total  = {total}")
    mean += (total / size).tolist()
    print(f"mean   = {mean}")
    numleft -= n
    print()
var = np.var(mean)
print(f"var    = {var}")

実行結果

sample
[[1]
 [0]
 [0]
 [0]
 [1]]
sum
[2]

sample
[[1]
 [0]
 [1]]
sum
[2]

total  = [4.]
mean   = [0.5]

sample
[[0]
 [0]
 [1]
 [0]
 [1]]
sum
[2]

sample
[[1]
 [0]
 [1]]
sum
[2]

total  = [4.]
mean   = [0.5, 0.5]

sample
[[1]
 [1]
 [1]
 [1]
 [0]]
sum
[4]

sample
[[0]
 [0]
 [1]]
sum
[1]

total  = [5.]
mean   = [0.5, 0.5, 0.625]

var    = 0.0034722222222222225

実行結果から3 つの標本が 1 つずつ順番に計算されて mean の要素に追加されていく様子が確認できます。

下記は maxe = 20 とした場合のプログラムです。この場合は maxs = min(8, 20) = 8 なので、一度の抽出で標本サイズが 8 の標本を抽出します。

また、maxn = min(3, 20 // 8) = min(3, 2) = 2 なので、一度の抽出で最大 2 つの標本を抽出します。実行結果から実際にそのような抽出が行われていることを確認して下さい。

size = 8
num = 3
maxe = 20

先程と同じなので略

var = np.var(mean)
print(f"var    = {var}")

省略していないプログラム

size = 8
num = 3
maxe = 20

numleft = num
maxs = min(size, maxe)
maxn = min(num, maxe // maxs)
mean = []
while numleft > 0:
    sizeleft = size
    n = min(numleft, maxn)
    total = np.zeros(n)
    while sizeleft > 0:
        s = min(sizeleft, maxe)
        sample = np.random.choice(X, p=P, size=(s, n)) 
        print("sample")
        print(f"{sample}")
        print("sum")
        print(f"{np.sum(sample, axis=0)}")
        print()
        total += np.sum(sample, axis=0)  
        sizeleft -= s
    print(f"total  = {total}")
    mean += (total / size).tolist()
    print(f"mean   = {mean}")
    numleft -= n
    print()
var = np.var(mean)
print(f"var    = {var}")

実行結果

sample
[[0 0]
 [1 0]
 [0 1]
 [1 0]
 [0 0]
 [1 0]
 [0 0]
 [0 1]]
sum
[3 2]

total  = [3. 2.]
mean   = [0.375, 0.25]

sample
[[1]
 [1]
 [1]
 [1]
 [0]
 [1]
 [0]
 [1]]
sum
[6]

total  = [6.]
mean   = [0.375, 0.25, 0.75]

var    = 0.045138888888888895

上記の 24 行目では mean に代入された list に対して += 演算子で list の結合を行いますが、+= の右に標本平均の一覧を表す total / size をそのまま記述すると、list と ndarray を + 演算子で計算するという意味になるため、list の結合とは異なる演算が行われるという問題が発生します。

ndarray にはその値を list に変換する tolist というメソッドがあるので、それを利用して ndarray を list に変換することで + 演算子による list の結合を行うことができるようになります。

なお、同様の処理を組み込み関数 list を利用した list(total / size) で list への変換を行うこともできますが、下記のプログラムのように組み込み関数 list では 2 次元以上の ndarray をうまく list に変換できないという点などが異なります。

a = np.array([[1, 2], [3, 4]])
b = list(a)
c = a.tolist()
print(b)
print(c)

実行結果

[array([1, 2]), array([3, 4])]
[[1, 2], [3, 4]]

numpy tolist メソッドの詳細については下記のリンク先を参照して下さい。

sampling_and_analyze の修正

下記は上記の分割処理を行うように sampling_and_analyze を修正したプログラムで、デフォルト値を 100 万とした maxe を代入する仮引数を追加しました。その他の修正点は標本抽出と標本平均の一覧の計算を上記で説明した方法で変更した点だけなので説明は省略します。ただし、確認のための print はもう必要がなくなったので削除しました。

def sampling_and_analyze(X, P, size, num, maxe=1000000):
    print(f"標本サイズ = {size}")
    E = np.average(X, weights=P)
    V = np.average(X ** 2, weights=P) - E ** 2
    print(f"母平均           = {E:.7f}")
    print(f"母分散           = {V:.7f}")

    numleft = num
    maxs = min(size, maxe)
    maxn = min(num, maxe // maxs)
    mean = []
    while numleft > 0:
        sizeleft = size
        n = min(numleft, maxn)
        total = np.zeros(n)
        while sizeleft > 0:
            s = min(sizeleft, maxe)
            sample = np.random.choice(X, p=P, size=(s, n)) 
            total += np.sum(sample, axis=0)  
            sizeleft -= s
        mean += (total / size).tolist()
        numleft -= n
    var = np.var(mean)
    print(f"標本平均の分散    = {var:.7f}") 
    print(f"母分散÷標本サイズ = {V/size:.7f}") 
    print(f"誤差             = {np.abs(var - V/size):.7f}") 

上記のプログラムでは標本平均の一覧を表す mean に標本数を表す num 個の要素が蓄積されていくので、num を例えば 1 億のような非常に大きな値にするとメモリが足りなくなる可能性が生じます。ただし、標本平均の分散を計算したり、この後で標本平均の分布が中心極限定理によって正規分布に近づくことを示す際に必要な標本数は数万もあれば十分なので、そのことは考慮しないことにします。

下記は sampling_and_analyze で下記の計算を行うプログラムです。

• X、P に偏りのないサイコロの出目を表す確率変数の値とその確率を代入する
• 標本サイズが 1, 10, 100, 1000, 10000 のそれぞれの場合に対して計算を行う
• 標本数を 10000 とする

実行結果から標本平均の分散が標本サイズに反比例し、いずれの場合も標本平均の分散と母分散÷標本サイズがほぼ等しくなることが確認できました。

X = np.array([1, 2, 3, 4, 5, 6])
P = np.array([1/6] * 6)
for size in [1, 10, 100, 1000, 10000]:
    sampling_and_analyze(X, P, size, 10000)
    print()

実行結果

標本サイズ = 1
母平均           = 3.5000000
母分散           = 2.9166667
標本平均の分散    = 2.9379211
母分散÷標本サイズ = 2.9166667
誤差             = 0.0212544

標本サイズ = 10
母平均           = 3.5000000
母分散           = 2.9166667
標本平均の分散    = 0.2912838
母分散÷標本サイズ = 0.2916667
誤差             = 0.0003828

標本サイズ = 100
母平均           = 3.5000000
母分散           = 2.9166667
標本平均の分散    = 0.0288969
母分散÷標本サイズ = 0.0291667
誤差             = 0.0002698

標本サイズ = 1000
母平均           = 3.5000000
母分散           = 2.9166667
標本平均の分散    = 0.0029491
母分散÷標本サイズ = 0.0029167
誤差             = 0.0000324

標本サイズ = 10000
母平均           = 3.5000000
母分散           = 2.9166667
標本平均の分散    = 0.0002873
母分散÷標本サイズ = 0.0002917
誤差             = 0.0000044

下記は以前の記事で定義した create_pd で作成した 1 から 10 までの整数の値をとる確率変数が従う確率分布に対して同様の計算を行うプログラムです。実行結果から、同様の結果になることが確認できました。

from util import create_pd

X, P = create_pd(minx=1, maxx=10)
print(X)
print(P)
for size in [1, 10, 100, 1000, 10000]:
    sampling_and_analyze(X, P, size, 10000)
    print()

実行結果

[ 1  2  3  4  5  6  7  8  9 10]
[0.03177847 0.16275253 0.13106408 0.11791219 0.15184416 0.04539966
 0.0634534  0.14413338 0.063114   0.08854814]
標本サイズ = 1
母平均           = 5.3044915
母分散           = 7.3426229
標本平均の分散    = 7.2980062
母分散÷標本サイズ = 7.3426229
誤差             = 0.0446166

標本サイズ = 10
母平均           = 5.3044915
母分散           = 7.3426229
標本平均の分散    = 0.7240390
母分散÷標本サイズ = 0.7342623
誤差             = 0.0102233

標本サイズ = 100
母平均           = 5.3044915
母分散           = 7.3426229
標本平均の分散    = 0.0745378
母分散÷標本サイズ = 0.0734262
誤差             = 0.0011116

標本サイズ = 1000
母平均           = 5.3044915
母分散           = 7.3426229
標本平均の分散    = 0.0073759
母分散÷標本サイズ = 0.0073426
誤差             = 0.0000333

標本サイズ = 10000
母平均           = 5.3044915
母分散           = 7.3426229
標本平均の分散    = 0.0007425
母分散÷標本サイズ = 0.0007343
誤差             = 0.0000082

中心極限定理の確認

次に、先程紹介した下記の中心極限定理が正しいことをグラフを描画することで確認することにします。

有限な値の母平均 $μ$ と母分散 $σ^2$ を持つ母集団の確率分布を表す確率変数を $X$ とし、その母集団から $n$ 個の要素を無作為復元抽出を行うものとする。作成された標本サイズが $n$ の標本の標本平均を表す確率変数を $\bar{X_n}$ とすると、$n$ を十分大きくすると $\bar{X_n}$ が従う確率分布は期待値が $μ$、分散が $\frac{σ ^ 2}{n}$ の正規分布で近似できる。

標本平均の分布を表す棒グラフの描画

標本平均の分布を表すグラフとして、それぞれの標本平均の値とその個数を表す棒グラフを描画するという方法を考えることができますが、そのようなグラフはあまり良いグラフにはなりません。そのことを最初に示します。

それぞれの標本平均の値とその個数は、前回の記事で説明したように np.unique で計算することができます。なお、np.bar で描画される棒グラフは、それぞれの棒の幅が初期設定では 1 となるため、それぞれの棒の幅を element の最大値と最小値の差を element の要素の個数で割った値にする必要があります。

厳密には特定の値がたまたま出現しなかった場合は element の要素の間隔が不均一になります。その結果、上記の方法で棒グラフの棒の幅を計算して描画すると、この下の標本サイズが 10 のグラフのように一部で棒グラフの隙間が空いたグラフが描画される場合があります。

なお、下記のプログラムは棒グラフによる可視化の限界（うまくいかない例）を体験してもらうためのプログラムなので、その点は気にせずに読み進めてください。

下記は先ほどの sampling_and_analyze を、標本平均のそれぞれの値の個数を表す棒グラフを描画するように修正したプログラムです。

• 5 行目：標本平均の一覧を表す element から np.unique を利用してそれぞれの標本の値とその数を表す element と count を計算する
• 6、7 行目：上記で説明した式で棒グラフの幅を計算し、棒グラフを描画する

1  import matplotlib.pyplot as plt
2  import japanize_matplotlib
3  
4  def sampling_and_analyze(X, P, size, num, maxe=1000000):
元と同じなので省略
5      element, count = np.unique(mean, return_counts=True) 
6      width = ((max(element) - min(element)) / len(element))
7      plt.bar(element, count, width=width)
8      plt.title(f"標本サイズ {size}")
9      plt.show()

行番号のないプログラム

import matplotlib.pyplot as plt

def sampling_and_analyze(X, P, size, num, maxe=1000000):
    print(f"標本サイズ = {size}")
    E = np.average(X, weights=P)
    V = np.average(X ** 2, weights=P) - E ** 2
    print(f"母平均           = {E:.7f}")
    print(f"母分散           = {V:.7f}")

    numleft = num
    maxs = min(size, maxe)
    maxn = min(num, maxe // maxs)
    mean = []
    while numleft > 0:
        sizeleft = size
        n = min(numleft, maxn)
        total = np.zeros(n)
        while sizeleft > 0:
            s = min(sizeleft, maxe)
            sample = np.random.choice(X, p=P, size=(s, n)) 
            total += np.sum(sample, axis=0)  
            sizeleft -= s
        mean += (total / size).tolist()
        numleft -= n
    var = np.var(mean)
    print(f"標本平均の分散    = {var:.7f}") 
    print(f"母分散÷標本サイズ = {V/size:.7f}") 
    print(f"誤差             = {np.abs(var - V/size):.7f}") 
    
    element, count = np.unique(mean, return_counts=True) 
    width = ((max(element) - min(element)) / len(element))
    plt.bar(element, count, width=width)
    plt.title(f"標本サイズ {size}")    
    plt.show()

下記は先ほどと同じ条件で sampling_and_analyze を呼び出すことで、標本平均を表す棒グラフを描画するプログラムです。母平均などの表示は先ほどと変わらないので、下記にはグラフだけを示します。

for size in [1, 10, 100, 1000, 10000]:
    sampling_and_analyze(X, P, size, 10000)
    print()

実行結果

　
　

標本平均は標本サイズが大きくなればなるほど取りうる値の種類が大きくなります。例えば 0 と 1 のいずれかの値を取る確率変数 $X$ から抽出された標本サイズが $n$ の標本の標本平均 $\bar{X_n}$ がとりうる値は $n = 1$ の場合は 0 と 1 の 2 種類ですが、$n = 2$ の場合は 0、0.5、1 の 3 種類、$n = 10$ の場合は 0、0.1、0.2、・・・、1 の 11 種類になります。確率変数 $X$ がとりうる値の種類が 2 より多い場合も同様です。

そのため、標本サイズが大きくなると標本平均がとる値の種類が大きく増えまてしまいます。一方で、sampling_and_analyze が抽出する標本数は変わらないので、細分化された個々の値の出現回数にランダムなばらつきが目立つようになり、上記の標本サイズが 1000 以上の場合の棒グラフのようにグラフがギザギザの形状になってしまいます。これが今回の検証で単純な棒グラフを利用しないほうが良い理由です。

標本平均の分布を表すヒストグラムの描画

標本の分布を表すグラフはそれぞれの標本の値の個数をグラフ化するのではなく、特定の範囲の個数をグラフ化するという、ヒストグラム（度数分布図）でグラフ化するのが一般的です。ヒストグラムでは個数¹を計算する範囲のことを区間または階級と呼び、英語では bin と呼びます。

先ほどの棒グラフはすべての実現値に対して個別に階級（ビン）を割り当てた、極端に階級数を多くしたヒストグラムと考えることができます。

matplotlib では hist という関数でヒストグラムを描画することができ、その際に bins という仮引数に階級の個数を代入します。下記は sampling_and_analyze にヒストグラムの階級の数を表す bins という仮引数を追加してヒストグラムを描画するように修正したプログラムです。

• 1 行目：ヒストグラムの階級の数を代入する仮引数 bins を追加した
• 3 行目：棒グラフの描画を行うプログラムを削除し、plt.hist で mean のデータを階級数が bins のヒストグラムで描画する

1  def sampling_and_analyze(X, P, size, num, bins, maxe=1000000):
元と同じなので省略
2  
3      plt.hist(mean, bins=bins)
4      plt.title(f"標本サイズ {size}")
5      plt.show()

行番号のないプログラム

def sampling_and_analyze(X, P, size, num, bins, maxe=1000000):
    print(f"標本サイズ = {size}")
    E = np.average(X, weights=P)
    V = np.average(X ** 2, weights=P) - E ** 2
    print(f"母平均           = {E:.7f}")
    print(f"母分散           = {V:.7f}")

    numleft = num
    maxs = min(size, maxe)
    maxn = min(num, maxe // maxs)
    mean = []
    while numleft > 0:
        sizeleft = size
        n = min(numleft, maxn)
        total = np.zeros(n)
        while sizeleft > 0:
            s = min(sizeleft, maxe)
            sample = np.random.choice(X, p=P, size=(s, n)) 
            total += np.sum(sample, axis=0)  
            sizeleft -= s
        mean += (total / size).tolist()
        numleft -= n
    var = np.var(mean)
    print(f"標本平均の分散    = {var:.7f}") 
    print(f"母分散÷標本サイズ = {V/size:.7f}") 
    print(f"誤差             = {np.abs(var - V/size):.7f}") 
    
    plt.hist(mean, bins=bins)
    plt.title(f"標本サイズ {size}")
    plt.show()

下記は sampling_and_analyze で標本サイズが 1000 の場合に、階級の数が 5、10、20、50、100 の場合のヒストグラムを描画するプログラムです。実行結果から階級の数によってヒストグラムの形状が変化し、見た目の印象が大きく変わることが確認できます。

for bins in [5, 10, 20, 50, 100]:
    sampling_and_analyze(X, P, 1000, 10000, bins=bins)

実行結果

　
　

参考までにヒストグラムの Wikipedia のリンクを下記に紹介します。

plt.hist の詳細は下記のリンク先を参照して下さい。

ヒストグラムの階級の数の調整

ヒストグラムの階級の数をいくつにすれば良いかについては統計学の世界でも様々な議論が行われていますが、どのようなデータに対しても常にふさわしい階級の数を計算する方法は残念ながら存在しません。本記事では試行錯誤の末、階級の数を下記のうちの最大値で計算することにしました。他にもっと良い階級の数を計算する方法を思いついた方はその方法でヒストグラムを描画してみて下さい。

• 確率分布がとりうる値の個数
• 「標本平均がとる値の種類の平方根 * 2」と 50 の最小値

階級の数の最小値を「確率分布がとりうる値の個数」とした理由は、標本サイズを 1 とした場合のヒストグラムの形状が確率分布を表すグラフの形状と一致するようにしたかったからです。なお、標本サイズを 1 とした場合の標本平均を表す確率変数 $\bar{X_1}$ はその定義から元の確率分布を表す確率変数 $X$ と等しくなります。

下記はそのように sampling_and_analyze を修正したプログラムです。

• 1 行目：仮引数 bins を None をデフォルト値とするデフォルト引数に修正した
• 3 ~ 6 行目：bins が None でない場合に上記で説明した値を計算して bins に代入するように修正した。なお、標本平均がとる値の種類は np.unique で計算された標本平均がとる値の一覧の個数を組み込み関数 len で計算することで求められる

1  def sampling_and_analyze(X, P, size, num, bins=None, maxe=1000000):
元と同じなので省略
2  
3      if bins is None:
4          element = np.unique(mean)
5          bins = max(len(X), min(int(np.sqrt(len(element))) * 2, 50))
6      plt.hist(mean, bins=bins)
7      plt.title(f"標本サイズ {size}")
8      plt.show()

行番号のないプログラム

def sampling_and_analyze(X, P, size, num, bins=None, maxe=1000000):
    print(f"標本サイズ = {size}")
    E = np.average(X, weights=P)
    V = np.average(X ** 2, weights=P) - E ** 2
    print(f"母平均           = {E:.7f}")
    print(f"母分散           = {V:.7f}")

    numleft = num
    maxs = min(size, maxe)
    maxn = min(numleft, maxe // maxs)
    mean = []
    while numleft > 0:
        sizeleft = size
        n = min(numleft, maxn)
        total = np.zeros(n)
        while sizeleft > 0:
            s = min(sizeleft, maxe)
            sample = np.random.choice(X, p=P, size=(s, n)) 
            total += np.sum(sample, axis=0)  
            sizeleft -= s
        mean += (total / size).tolist()
        numleft -= n
    var = np.var(mean)
    print(f"標本平均の分散    = {var:.7f}") 
    print(f"母分散÷標本サイズ = {V/size:.7f}") 
    print(f"誤差             = {np.abs(var - V/size):.7f}") 
    
    if bins is None:
        element = np.unique(mean)
        bins = max(len(X), min(int(np.sqrt(len(element))) * 2, 50))
    plt.hist(mean, bins=bins)
    plt.title(f"標本サイズ {size}")
    plt.show()

下記は修正した sampling_and_analyze でヒストグラムを描画するプログラムです。なお、標本数が 10000 の場合はばらつきが多いためあまりきれいなヒストグラムが描画されないことが分かったので標本数を 100000 としました。実行結果から、見た目がきれいなヒストグラムが描画できているのではないかと思います。

for size in [1, 10, 100, 1000, 10000]:
    sampling_and_analyze(X, P, size, 100000)
    print()

実行結果

　
　

参考までにヒストグラムの階級の幅と個数に関する Wikipedia の項目のリンクを下記に示します。

正規分布を表すグラフの描画

正規分布を表す確率変数は連続型確率変数で、平均が m、標準偏差が s である正規分布の確率分布を表す確率密度関数は以前の記事で説明したように statistics モジュールの NormalDist というクラスを利用して計算することができます。

確率密度関数は、そのグラフと x 軸の間の面積が 1 になるという性質があります。下図は平均が 0、標準偏差が 1 の正規分布のグラフを描画するプログラムで、赤色の部分の面積が 1 になります。なお、このプログラムの意味については以前の記事を参照して下さい。

from statistics import NormalDist

ndist = NormalDist(0, 1)
NX = [x / 100 for x in range(-500, 500)]
NY = [ndist.pdf(x) for x in NX]
plt.plot(NX, NY)
plt.fill_between(NX, 0, NY, color="r")
plt.show()

実行結果

また、上記の正規分布を表すグラフの縦軸の数値は 1 未満の小さな値になります。そのため、上記の正規分布を表すグラフと 1 よりはるかに大きな階級の個数をグラフ化するヒストグラムを重ねて描画しても縦軸の縮尺（スケール）が違いすぎてうまく比較できないという問題があります。

この問題を解決するために、ヒストグラムを描画する plt.hist には、確率密度関数と同様にグラフの面積の合計が 1 になるように描画を行うこと指定する仮引数 density² が用意されています。density に True を代入して描画することでヒストグラムの縦軸の縮尺が確率密度関数と一致し、2 つのグラフを重ねて比較できるようになります。

そこで、sampling_and_analyze を下記のように修正することにします。

• plt.hist をキーワード引数 density=True を記述して呼び出すようにする
• ヒストグラムを描画する範囲にあわせて中心極限定理の理論上の曲線を表す、平均が $μ$、分散が $\frac{σ^2}{n}$ の正規分布のグラフを重ねて描画する

下記は sampling_and_analyze をそのように修正したプログラムです。

• 8 行目：平均が E、標準偏差が np.sqrt(V/size) の正規分布を表す NormalDist クラスのインスタンスを作成する。なお、NormalDist では正規分布の標準偏差を指定する必要があり、標準偏差は分散の平方根で計算できる
• 9 行目：以前の記事で紹介した、np.linspace で標本平均の最小値と最大値の間を 100 分割した値を要素とする ndarray を計算して NX に代入する
• 10 行目：ndist.pdf を利用して NXの各値に対する確率密度関数の値を計算し、それらを要素とする listを NY に代入する
• 11 行目：plt.plot の実引数に NX、NY を記述して理論上の正規分布を表す曲線のグラフを重ねて描画する

 1  def sampling_and_analyze(X, P, size, num, bins=None, maxe=1000000):
元と同じなので省略 
 2  
 3      if bins is None:
 4          element = np.unique(mean)
 5          print(len(element))
 6          bins = max(len(X), min(int(np.sqrt(len(element))) * 2, 50))
 7      plt.hist(mean, bins=bins, density=True)
 8      ndist = NormalDist(E, np.sqrt(V/size))
 9      NX = np.linspace(min(mean), max(mean), 50)
10      NY = [ndist.pdf(x) for x in NX]
11      plt.plot(NX, NY)
12      plt.title(f"標本サイズ {size}")
13      plt.show()

行番号のないプログラム

def sampling_and_analyze(X, P, size, num, bins=None, maxe=1000000):
    print(f"標本サイズ = {size}")
    E = np.average(X, weights=P)
    V = np.average(X ** 2, weights=P) - E ** 2
    print(f"母平均           = {E:.7f}")
    print(f"母分散           = {V:.7f}")

    numleft = num
    maxs = min(size, maxe)
    maxn = min(numleft, maxe // maxs)
    mean = []
    while numleft > 0:
        sizeleft = size
        n = min(numleft, maxn)
        total = np.zeros(n)
        while sizeleft > 0:
            s = min(sizeleft, maxe)
            sample = np.random.choice(X, p=P, size=(s, n)) 
            total += np.sum(sample, axis=0)  
            sizeleft -= s
        mean += (total / size).tolist()
        numleft -= n
    var = np.var(mean)
    print(f"標本平均の分散    = {var:.7f}") 
    print(f"母分散÷標本サイズ = {V/size:.7f}") 
    print(f"誤差             = {np.abs(var - V/size):.7f}") 
    
    if bins is None:
        element = np.unique(mean)
        bins = max(len(X), min(int(np.sqrt(len(element))) * 2, 50))
    plt.hist(mean, bins=bins, density=True)
    ndist = NormalDist(E, np.sqrt(V/size))
    NX = np.linspace(min(mean), max(mean), 100)
    NY = [ndist.pdf(x) for x in NX]
    plt.plot(NX, NY)
    plt.title(f"標本サイズ {size}")
    plt.show()

下記は、これまでの条件に標本サイズが 2、3、5 の場合を加えて上記で修正した sampling_and_analyze を実行するプログラムです。実行結果から標本サイズが大きくなるにつれて、標本平均のヒストグラムの形状が理論上の正規分布の曲線に急速にきれいに重なっていく（近似できるようになる）様子がわかります。また、一般的に言われている標本サイズが 30 以上になる前のサイズが 5 や 10 の段階でも、かなり正規分布に近い形状になることが視覚的に実感できるはずです。

for size in [1, 2, 3, 5, 10, 100, 1000, 10000]:
    sampling_and_analyze(X, P, size, 100000)

実行結果

　
　
　
　

なお、標本サイズが一定以上になると描画されるグラフの形状がほぼ同一になりますが、実際にはそれらのグラフの横軸の縮尺（スケール）が大きく異なります。そのことは、標本サイズが大きくなるにつれてグラフの x 軸の範囲がどんどん狭くなっていること（大数の法則の性質）からも確認することができます。

上記から、中心極限定理が主張する「分布収束」の現象が本当に成り立っていることが視覚的に確認できました。興味がある方は create_pd を利用してさらに極端に偏った形状の確率分布などを作成し、どのような確率分布であっても最終的に標本平均の分布が正規分布に収束するという、中心極限定理の凄さをぜひ体験してみて下さい。

中心極限定理と区間推定

中心極限定理が非常に便利な点は、有限な期待値と分散を持つという条件を満たせば元の確率分布がどのような分布であっても、標本平均の確率分布を正規分布で近似できるということです。この性質を利用することで、無作為復元抽出を行って実際に得られた標本から、母集団の母平均がどの範囲にどれくらいの確率で含まれるかを推定するという「区間推定」を行うことができます。

正規分布の重要な性質として、平均を中心とした標準偏差の特定の倍数の間の範囲に、データがどれくらいの確率で収まるかを計算できるというものがあります。例えば、良く利用される基準として標準偏差の 1.96 倍、すなわち「平均 ± 1.96×標準偏差」の範囲にデータが収まる確率は 95 % になります。下記はそのことを図示するプログラムで、平均が 0、標準偏差が 1 の正規分布の $-1.96 ≦ x ≦ 1.96$ の範囲を表す黄色の部分の面積が全体の 95 % の 0.95 に、それ以外の赤色の部分の面積が残りの 5 % の 0.05 になります。

ndist = NormalDist(0, 1)
NX = [x / 100 for x in range(-500, 500)]
NY = [ndist.pdf(x) for x in NX]
plt.plot(NX, NY)
plt.fill_between(NX, 0, NY, color="r")
plt.fill_between(NX[500-196:500+196], 0, NY[500-196:500+196], color="y")
plt.show()

実行結果

中心極限定理により、標本サイズが $n$ の標本平均が従う分布の平均は $μ$、標準偏差は $\frac{σ}{\sqrt{n}}$ となります。従って、実際に得られた標本の標本平均を $\bar{x_n}$ とすると、95% の確率で下記の不等式が成り立ちます。

$μ - 1.96\frac{σ}{\sqrt{n}} ≦ \bar{x_n} ≦ μ + 1.96\frac{σ}{\sqrt{n}}$

下記はこの式を未知の母平均 $μ$ が中心となるように変形した式です。

$\bar{x_n} - 1.96\frac{σ}{\sqrt{n}} ≦ μ ≦ \bar{x_n} + 1.96\frac{σ}{\sqrt{n}}$

この式は、標本抽出によって実際に得られた標本から計算された上記の左右の式の範囲が 95% の確率で母平均 $μ$ を含んでいることを意味します。このような範囲のことを信頼区間と呼び、上記を「母平均の 95% 信頼区間」と呼びます。また、このように実際に得られた標本から母集団の統計量の信頼区間を推定することを区間推定と呼びます。

なお、上記の信頼区間の左右の式の中の $\bar{x_n}$ は実際に得られた標本から計算できますが、母分散 $σ^2$（および母分散から計算される母標準偏差 $σ$）は理論上の値なので、現実のデータ分析では事前に知ることはできません。そのため、母平均の区間推定では母分散を以前の記事で説明した標本から計算できる不偏分散で代用します。

誤差を含む不偏分散は母分散とは異なる値になるため、不偏分散で代用を行うと正規分布ではない t 分布と呼ばれる別の確率分布をベースにする必要があり、信頼区間の計算式はさらに複雑になります。ただし、区間推定を行う際にはこれらの複雑な数式を手作業で計算する必要はなく、統計学に関する処理が実装されたモジュールに用意されている区間推定の関数を利用して簡単に算出することができます。

区間推定の詳しいメカニズムや t 分布の性質については難しいため本記事では省略します。興味がある方は、下記の Wikipedia の項目などを参照して調べてみてください。

有限な値の分散を持たない確率分布に対する大数の法則

以下の内容はおまけです。

以前の記事で行った大数の弱法則の証明では確率分布が有限な値の分散を持つことを前提にしていましたが、大数の弱法則は確率分布が有限な値の分散を持たない場合でも成り立ちます。その証明は難しいので省略しますが、実際にそのことが成り立つことを具体例を挙げて示すことにします。

有限な値の分散を持たない確率分布の具体例

有限な期待値を持つが有限な分散を持たない確率変数は、以前の記事で説明したサンクトペテルブルクのパラドックスに出てくる期待値が有限な値を持たない確率変数を $X$ をベースにして、$Y = \sqrt{X}$ という式で定義することができるので最初にそのことを証明します。

下記はサンクトペテルブルクのパラドックスに出てくる期待値が有限な値を持たない確率分布 $X$ の確率質量関数です。

• $P(X=2^{n-1}) = (\frac{1}{2})^{n}$ （$n$ が自然数の場合）
• $P(X=x) = 0$ （$x$ が $2^{n-1}$ 以外の場合）

ここで新しい確率変数 $Y$ を $Y = \sqrt{X}$ という式で定義すると、確率変数 $Y$ の確率質量関数は下記のようになります。

• $P(Y=\sqrt{2^{n-1}}) = (\frac{1}{2})^{n}$ （$n$ 自然数の場合）
• $P(Y=y) = 0$（$y$ が $\sqrt{2^{n-1}}$ 以外の場合）

従って、確率変数 $Y$ の期待値 $E[Y]$ は下記の式で計算することができます。

$E[Y] = \sum_{n=1}^{∞} (\sqrt{2^{n-1}})(\frac{1}{2})^{n}$

$= \sum_{n=1}^{∞} (\sqrt{2})^{n-1}2^{-n}$

$= \sum_{n=1}^{∞} (\sqrt{2})^{n-1}(\sqrt{2}^2)^{-n}$

$= \sum_{n=1}^{∞} (\sqrt{2})^{n-1}(\sqrt{2})^{-2n}$

$= \sum_{n=1}^{∞} (\sqrt{2})^{-n-1}$

$= \sum_{n=1}^{∞} (\sqrt{2})^{-(n+1)}$

$= \sum_{n=1}^{∞} (\frac{1}{\sqrt{2}})^{n + 1}$

ここで $m = n - 1$ おくと $n = 1$ のときに $m = 0$、$n + 1 = m + 2$ となるので上記の式は下記のようになります。

$= \sum_{m=0}^{∞} (\frac{1}{\sqrt{2}})^{m+2}$

$= \sum_{m=0}^{∞} (\frac{1}{\sqrt{2}})^m(\frac{1}{\sqrt{2}})^2$

$= \sum_{m=0}^{∞} \frac{1}{2}(\frac{1}{\sqrt{2}})^m$

$\sqrt{2} = 約 1.41$ から $0 < \frac{1}{\sqrt{2}} < 1$ が成り立つので、以前の記事 で説明した無限等比数列の和の公式 $\sum_{n=0}^{∞} ar^i = \frac{a}{1-r}$ （ただし $-1 < r < 1$）を適用すると上記の式は下記のように約 1.707 という有限の値に収束します。

$= \frac{1}{2}\frac{1}{1-\frac{1}{\sqrt{2}}}$

$= \frac{1}{2}\frac{\sqrt{2}}{\sqrt{2} - 1}$

$= \frac{1}{2}\frac{\sqrt{2}(\sqrt{2} + 1)}{(\sqrt{2} - 1)(\sqrt{2} + 1)}$

$= \frac{1}{2}\frac{2 + \sqrt{2}}{2 - 1}$

$= 1 + \frac{\sqrt{2}}{2} = 約 1.707$

一方、確率変数 $Y$ の分散は以前の記事で説明したように下記の式で計算できます。

$V[Y] = E[Y^2] - E[Y]^2$

$Y = \sqrt{X}$ から $Y^2 = X$ となり、確率変数 $X$ が有限の期待値を持たないサンクトペテルブルクのパラドックスの確率分布に従うことから、確率変数 $Y^2$ の期待値である $E[Y^2] = E[X]$ は有限な値を持たないことがわかります。先ほど示したように $E[Y]$ は約 1.7 という有限な値ですが、$E[Y^2]$ は有限な値ではないので上記の式から $V[Y]$ は有限な値になりません。

上記から、確率変数 $Y$ は有限な値の期待値を持つが、有限な値の分散を持たないことが証明されました。

betY の定義

上記の確率変数 $Y$ に従う確率変数から無作為復元抽出を行い、大数の法則が成り立つことを検証する関数は、以前の記事で定義した確率変数 $X$ から無作為復元抽出を行う bet3 を改良することで簡単に定義することができます。

具体的には bet3 で 2 のべき乗を計算していた部分を、下記のプログラムのように $\sqrt{2}$（2 ** 0.5)のべき乗を計算するように修正するだけです。なお、下記ではこの関数の名前を betY としました。

def betY(size):
    prize = np.pow(2 ** 0.5, -np.floor(np.log2(np.random.uniform(size=size))) - 1)
    mean = np.cumsum(prize) / np.arange(1, size + 1)
    plt.plot(mean)
    plt.xlabel("賭けの回数")
    plt.ylabel("賞金の平均")   

下記は betY で 1000 万回無作為復元抽出を行い、試行回数を増やした場合の賞金³の累積平均の変化を折れ線グラフで描画するプログラムです。実行結果から以前の記事 で紹介したサンクトペテルブルクのパラドックスの確率分布と異なり、試行回数が多くなるにつれて累積平均の値が先ほど数学的に計算した理論上の期待値（平均）である 約 1.7 へと綺麗に収束していく様子が確認できます。

betY(10000000)

実行結果

今回の記事のまとめ

今回の記事では標本平均の分散の性質と中心極限定理が正しいことを視覚的に確認し、区間推定という中心極限定理の利用例について紹介しました。

また、おまけとして有限な期待値を持つが有限な分散を持たない確率分布が大数の法則を満たすことを視覚的に確認しました。

本記事で入力したプログラム

今回の記事で定義した関数は util.py に保存することにします。

リンク	説明
marubatsu.ipynb	本記事で入力して実行した JupyterLab のファイル
util.py	本記事で更新した util_new.py

次回の記事

近日公開予定です

1. 統計学では度数と呼びます。度数は度数分布表の名前の中に入っています ↩

2. density は密度を表す英単語です ↩

3. $n$ 回目で当たった場合の賞金額はサンクトペテルブルクのパラドックスの場合の $2^{n-1}$ とは異なり $\sqrt{2^{n-1}}$ です ↩