種々の確率分布関数 †
二項分布 †
% matplotlib inline
import math
import numpy as np
import matplotlib.pyplot as plt
def binomial(n,m,p):
comb = math.factorial(m) / (math.factorial(n) * math.factorial(m-n))
prob = comb * p ** n * (1 - p) ** (m - n)
return prob
p = 0.5 # １回の試行でPかQのどちらかが必ず起こる場合にPが起こる確率
M = 100 # 試行回数
N = 100000 # サンプル数
np.random.seed(0) # 乱数の初期化
R=np.random.binomial(n=M, p=p, size=N) # 乱数を用いてM回試行し，Pが起こった回数をRに保存（それをN回繰り返す）
plt.hist(R, bins=20,normed=True) # Rに格納された回数の分布を100本の棒グラフにする．さらに総面積が1になるよう規格化
x = np.arange(M) # 0からMの範囲を間隔1で刻み，配列xに保存
y = np.zeros(M) # 配列yの全要素を0に初期化
for i in range(M):
y[i]=binomial(i,M,p) # 二項分布の理論値を配列yに保存
plt.plot(x, y, lw=2, color='g') # x vs. yを線幅2の赤い線で描く
ax = plt.axes()
plt.title(r'$\mathrm{Binomial Distribution}\ p=%f,\ M=%i$' % (p,M))
plt.xlabel('$n$', fontsize=16)
plt.ylabel('$P(n;M)$', fontsize=16)
plt.show() # 上記のグラフをまとめて1つの図に出力
正規分布（Gauss分布） †
N = 100000 # 発生する乱数の数
std = 1.0
ave = 0.0
np.random.seed(0) # 乱数の初期化
R = std*np.random.randn(N)+ave # 平均=ave，標準偏差=stdの正規分布に従う乱数を発生させ，Rに格納
plt.hist(R, bins=100,normed=True) # Rに格納された乱数の分布を100本の棒グラフにする．さらに総面積が1になるよう規格化
x = np.arange(-5, 5, 0.01) # -5から5の範囲を間隔0.01で刻み，配列xに保存
y = np.exp(-(x-ave)**2/2/std**2)/np.sqrt(2*np.pi*std**2) # 配列xと同じ数だけ正規分布の理論値を配列yに保存
plt.plot(x, y, lw=2, color='r') # x vs. yを線幅2の赤い線で描く
plt.show() # 上記のグラフをまとめて1つの図に出力
中心極限定理の一例（二項分布→正規分布） †
def binomial(n,m,p):
comb = math.factorial(m) / (math.factorial(n) * math.factorial(m-n))
prob = comb * p ** n * (1 - p) ** (m - n)
return prob
p = 0.5 # １回の試行で事象が起こる確率
M = 100 # 試行回数
N = 100000 # サンプル数
np.random.seed(0) # 乱数の初期化
R=np.random.binomial(n=M, p=p, size=N) # 乱数を用いてM回試行し，pが起こった回数をRに保存（それをN回繰り返す）
plt.hist(R, bins=20,normed=True) # Rに格納された回数の分布を100本の棒グラフにする．さらに総面積が1になるよう規格化
x = np.arange(M) # 0からMの範囲を間隔1で刻み，配列xに保存
y = np.zeros(M) # 配列yの全要素を0に初期化
for i in range(M):
y[i]=binomial(i,M,p) # 二項分布の理論値を配列yに保存
plt.plot(x, y, lw=2, color='g') # x vs. yを線幅2の赤い線で描く
ax = plt.axes()
plt.title(r'$\mathrm{Binomial Distribution}\ p=%f,\ M=%i$' % (p,M))
plt.xlabel('$n$', fontsize=16)
plt.ylabel('$P(n;M)$', fontsize=16)
std = np.sqrt(M*p*(1-p)) # 正規分布の標準偏差を二項分布に合うように設定
ave = p*M # 正規分布の平均値を二項分布に合うように設定
y = np.exp(-(x-ave)**2/2/std**2)/np.sqrt(2*np.pi*std**2) # 配列xと同じ数だけ正規分布の理論値を配列yに保存
plt.plot(x, y, lw=2, color='r') # x vs. yを線幅2の赤い線で描く
plt.show() # 上記のグラフをまとめて1つの図に出力
ランダムウォーク（酔歩） †
M = 10000 # 酔歩の回数
np.random.seed(0) # 乱数の初期化
R = np.random.choice([-1,1],M) # +1 or -1 をM個生成し，R(i)に保存
plt.hist(R, bins=20,normed=True) # Rに格納された乱数の分布を100本の棒グラフにする．さらに総面積が1になるよう規格化
plt.show() # 上記のグラフを出力
M = 10000 # 発生する乱数の数
L = 10000 # 独立なサンプル数
R = np.zeros(L)
np.random.seed(0) # 乱数の初期化
for i in range(L): # L回繰り返す
step = np.random.choice([-1,1],M) # +1/-1 をN個生成
R[i] = np.sum(step) # i回目のサンプルについて，N個発生させた+1/-1の乱数の合計値のをR(i)に保存
plt.hist(R, bins=20,normed=True) # Rに格納された乱数の分布を20本の棒グラフにする．さらに総面積が1になるよう規格化
p=0.5
std = 2*np.sqrt(M*p*(1-p)) # 正規分布の標準偏差を二項分布に合うように設定
ave = 0. # 正規分布の平均値を二項分布に合うように設定
x = np.arange(-400,400,10) # 配列xを作成
y = np.exp(-(x-ave)**2/2/std**2)/np.sqrt(2*np.pi*std**2) # 配列xと同じ数だけ正規分布の理論値を配列yに保存
plt.plot(x, y, lw=2, color='r') # x vs. yを線幅2の赤い線で描く
plt.show() # 上記のグラフを出力
宿題 †
- ランダムウォーク(+1/-1)をM回繰り返した後の位置iの分布関数（本来は二項分布P((M+i)/2;M)に従う）は，Mが大きい時には，平均<i>=0，分散σ^2=4xMx0.5^2の正規分布で近似することができる（M=10^4については確認済み）．M=10^2, 10^3, 10^4, 10^5の場合について，これを確かめよ．
