数学の具体的な問題にPythonを使って、数学もPythonも同時に学んでしまいましょう。今回はPythonを使った確率・統計の問題として、シューハートのノーマル・チップスをみてみたいと思います。シューハートのノーマル・チップスは、0から60までの数字の書かれたチップを998枚用意して1枚のチップをランダムに取り出すとき、そのチップに書かれた数字によって近似的に正規分布にしたがう確率変数を実現するものです。
シューハートのノーマル・チップス
0から60までの数字の書かれたカード(チップ)を998枚を以下のような枚数で準備します。
数字 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
枚数 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 3 | 4 | 4 | 5 | 7 | 8 | 9 | 11 | |
15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 | |
13 | 15 | 17 | 19 | 22 | 24 | 27 | 29 | 31 | 33 | 35 | 37 | 38 | 39 | 40 | 40 | |
31 | 32 | 33 | 34 | 35 | 36 | 37 | 38 | 39 | 40 | 41 | 42 | 43 | 44 | 45 | ||
40 | 39 | 38 | 37 | 35 | 33 | 31 | 29 | 27 | 24 | 22 | 19 | 17 | 15 | 13 | ||
46 | 47 | 48 | 49 | 50 | 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | ||
11 | 9 | 8 | 7 | 5 | 4 | 4 | 3 | 2 | 2 | 1 | 1 | 1 | 1 | 1 |
各数字の分布状況は以下の散布図で視覚的に捉えることができます。
chips = list(range(61)) chips_wheight = [1,1,1,1,1,2,2,3,4,4,5,7,8,9,11, 13,15,17,19,22,24,27,29,31,33,35,37,38,39,40,40, 40,39,38,37,35,33,31,29,27,24,22,19,17,15,13, 11,9,8,7,5,4,4,3,2,2,1,1,1,1,1] import matplotlib.pyplot as plt plt.scatter(chips, chips_wheight) plt.show()
より
という分布を得ます。一見して、平均値が30、標準偏差10の正規分布に近い形をしていることがわかります。以下、実際に平均値と標準偏差を計算して、これを確認してみます。
平均値・標準偏差
まず平均値(expectation value)と分散(variance)を定義通りに計算しましょう。まず平均値は
expec = [x*y for (x, y) in zip(chips, chips_wheight)] sum(expec)/sum(chips_wheight)
より
がわかります。次に分散ですが
varian = [(x-30.0)**2*y for (x, y) in zip(chips, chips_wheight)] sum(varian)/sum(chips_wheight)
より
がわかります。よって標準偏差は
import math math.sqrt(99.09)
より
を得ます。
正規分布
平均値 標準偏差 の正規分布は次の確率密度関数で与えられます。
関数をプロットしてみましょう。
import numpy as np mu = 30.0 sigma = 9.95 x = np.linspace( -10, 70, 400) y = 1/np.sqrt(2*np.pi*sigma**2)*np.exp( -(x-mu)**2/(2*sigma**2) ) plt.plot(x, y) plt.show()
より
という分布を得ます。
正規分布との比較
ノーマル・チップスの分布図との比較のために、正規分布をチップの枚数998倍だけして、合わせてプロットしてみましょう。
import numpy as np mu = 30.0 sigma = 9.95 x = np.linspace( -10, 70, 400) y = 998*1/np.sqrt(2*np.pi*sigma**2)*np.exp( -(x-mu)**2/(2*sigma**2) ) plt.plot(x, y) plt.scatter(chips, chips_wheight) plt.show()
より
となって、ノーマル・チップスが正規分布をよく近似していることがわかります。
今回は、シューハートのノーマル・チップスを解説してみました。現在では正規分布に従う確率変数はコンピュータで簡単に実現できますが、このような比較的シンプルな設定で正規分布が近似的に実現できるのは面白いですね。本記事の執筆者がシューハートのノーマル・チップスを知ったのは、以下に紹介している『統計学のすすめ 決定と計画の科学』(佐藤信・ブルーバックス)の中でした。古い本ですが、Kindleでも入手でき、数式の詳細に深入りせずに統計学の基本的な考え方を解説している良書です。