Pythonプログラミング(ステップ7・乱数と中心極限定理)
このページでは、確率や統計の分野で極めて重要な中心極限定理について考える。
1. 擬似乱数
コンピュータを使って、サイコロ投げなどのランダムな事象を模倣する際に、乱数を自動発生させる機能があると便利だ。 さもないと、「人間がサイコロを振って、出た目を、パソコンに入力する」といった作業を強いられることになる。 そのため、Pythonには乱数を発生させるための関数が用意されている。
使い方は、まず宣言部に
import random
と書いてから、random.randint(n,m)
と呼ぶと、nからm(m自身を含む)までの間の整数が
等しい確率で、ランダムに返される。つまり、呼び出す度に、値が変わる。
また、random.random()
と呼ぶと、0から1の範囲(1は含まない)のランダムな実数が返される。
擬似乱数を使って、「サイコロ」を10万回振り、出た目(リストの添字との対応がつけやすいので、サイコロの目は0から5まで、と考えた)の数と割合を出力するプログラム例を以下に示す:
右のコードでサイコロを振るところは
i=random.choice([0,1,2,3,4,5])
と書くこともできる。
# coding: utf-8 import random dice=[0,0,0,0,0,0] for cnt in range(0,100000,1): i = random.randint(0,5) dice[i] = dice[i] + 1 for i in range(0,6,1): print("サイコロの目:",i," 回数:",dice[i]," 確率:",dice[i]/cnt)
実数値の乱数
物理シミュレーション等では、特定の範囲の実数の乱数が必要になる場合が多い。aからbまでの範囲の実数乱数を得るには、0から1までの一様乱数を使って
x = a + random.random() * (b-a)
のように記述しても良いし、
x = random.uniform(a,b)
でもよい。ただし、前者は上端bを含まないのに対して、後者はbを含む点で違いがある。
練習:モンテカルロ法
乱数を用いた算法は、一般にモンテカルロ法と呼ばれている。
以下の考え方に基づいて、シミュレーションによって円周率の値を推定するプログラムを作成してみなさい。
1: 「大きな」数 n を決める 2: カウンター c をリセットする 3: 5までをn回反復する: 4: 0から1の間で一様に分布する二つの乱数 x, y を生成する 5: もし x*x + y*y が1より小さければ((x,y)が半径1の円の内部ならば)カウンター c を 1 だけ増やす 6: nが十分大きければ、この時点で、c/nは半径1の円の面積の1/4 (すなわちπ/4)になっているはずである。
ヒント
サンプル数を
2. 乱数の和
では、次に、
そこで、実際に、コンピュータ上でサイコロを振って、この様子を確かめてみた。
下のグラフは、『サイコロを
すると、よく見慣れたベル型の分布が現れるのが分かる。
この様子を数学的に表現したのが中心極限定理である。
すなわち、
上のグラフからも分かるように、元の
上記の性質を使うと、一様乱数
サイコロを3,5,10回振って目の合計を調べる試行を10000回行い、頻度分布をまとめてプロットするPythonコードの例。
# coding: utf-8 import random import numpy as np import matplotlib.pyplot as plt ntrial = 10000 dice3=[ ] dice5=[ ] dice10=[ ] for _ in range(ntrial): s3=0 for _ in range(3): s3 = s3 + random.randint(1,6) dice3.append(s3) s5=0 for _ in range(5): s5 = s5 + random.randint(1,6) dice5.append(s5) s10=0 for _ in range(10): s10 = s10 + random.randint(1,6) dice10.append(s10) bins = np.linspace(0, 100, 101) plt.hist([dice3,dice5,dice10], bins, label=["n=3","n=5","n=10"]) plt.xlabel('sum') plt.ylabel('Frequency') plt.legend(loc='upper right') plt.xlim(0,60) plt.grid(True) plt.show()
練習:正規乱数の生成
上記のアルゴリズムに従い、1〜6の目を持つサイコロの値を加えることで、正規分布
ヒント
サイコロの場合、平均