ry/Ilas2_04 - Soft Matter Engineering Gr.

種々の確率分布関数 †

必要なライブラリの読み込み †

% matplotlib inline
import math
import numpy as np
import matplotlib.pyplot as plt

↑

正規分布（Gauss分布） †

"np.random.randn(N)"を用いて，平均&mathjax{\langle R\rangle=0};，分散&mathjax{\sigma^2=1};の正規分布に従う乱数&mathjax{R_i};をN個発生する．

N = 100000                           # 発生する乱数の数
std = 1.0                            # 発生する乱数の標準偏差
ave = 0.0                            # 発生する乱数の平均値
np.random.seed(0)                    # 乱数の初期化
R = std*np.random.randn(N)+ave       # 平均=ave，標準偏差=stdの正規分布に従う乱数をN個発生させ，Rに格納
plt.xlabel('$i$', fontsize=16)
plt.ylabel('$R_i$', fontsize=16)
plt.xlim(0,1000)
plt.ylim(-5,5)
plt.plot(R)                          # R_i vs. i (i=1,2,...,N)
plt.show()

乱数&mathjax{R_i};の分布が正規分布&mathjax{P(R)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left[-\frac{(R-\langle R \rangle)^2}{2\sigma^2}\right]};に一致することを確認する．

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.xlabel('$R$', fontsize=16)
plt.ylabel('$P(R)$', fontsize=16)
plt.legend(['Gauss','histgram'], fontsize=14)
plt.show()                           # 上記のグラフをまとめて1つの図に出力

乱数&mathjax{R_i};の自己相関関数&mathjax{C(k)=\frac{1}{N}\sum_{i=1}^{N} (R_i - \langle R\rangle)(R_{i+k}-\langle R\rangle)};を計算し，乱数間に相関がないことを確認する．

def auto_correlate(dt):
    cor = np.correlate(dt,dt,mode="full")
    return cor[N-1:]
    return cor
corr = np.zeros(N)
corr=auto_correlate(R-ave)/N
plt.plot(corr)
plt.xlim(-100,1000)
plt.ylim(-0.2,1.2)
plt.xlabel('$k$', fontsize=16)
plt.ylabel('$C(k)$', fontsize=16)
plt.show()

↑

二項分布 †

"np.random.binomial(n=M, p=p, size=N)"を用いて，二項分布&mathjax{P(R,M)};に従う乱数&mathjax{R_i};をN個発生する．

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.xlabel('$i$', fontsize=16)
plt.ylabel('$R_i$', fontsize=16)
plt.xlim(0,1000)
#plt.ylim(-5,5)
plt.plot(R)
plt.show()

乱数&mathjax{R_i};の分布が二項分布&mathjax{P(R;M)=\frac{M!}{R!(M-R)!}p^{R}(1-p)^{M-R}};に一致することを確認する．

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'Binomial Distribution: p=%f, M=%i' % (p,M))
plt.xlabel('$R$', fontsize=16)
plt.ylabel('$P(R;M)$', fontsize=16)
ax.legend(['binomial','histgram'], fontsize=14)
plt.show()                                             # 上記のグラフをまとめて1つの図に出力

乱数&mathjax{R_i};の自己相関関数&mathjax{C(k)=\frac{1}{N}\sum_{i=1}^{N} (R_i - \langle R\rangle)(R_{i+k}-\langle R\rangle)};を計算し，乱数間に相関がないことを確認する．

def auto_correlate(dt):
    cor = np.correlate(dt,dt,mode="full")
    return cor[N-1:]
    return cor
corr = np.zeros(N)
ave = M*p
corr=auto_correlate(R-ave)/N
plt.plot(corr)
plt.xlim(-100,1000)
#plt.ylim(-0.2,1.2)
plt.xlabel('$k$', fontsize=16)
plt.ylabel('$C(k)$', fontsize=16)
plt.show()

↑

中心極限定理の一例（二項分布→正規分布） †

試行回数&mathjax{M};と確率&mathjax{p};が十分に大きいとき，二項分布&mathjax{P(R;M)};は平均&mathjax{\langle R\rangle=Mp};，分散&mathjax{\sigma^2=Mp(1-p)};の正規分布で近似できることを確認する．

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=3, color='g')         # x vs. yを線幅2の赤い線で描く
ax = plt.axes()
plt.title(r'Binomial Distribution: p=%f, M=%i' % (p,M))
plt.xlabel('$R$', fontsize=16)
plt.ylabel('$P(R;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の赤い線で描く
ax.legend(['binomial','Gauss','histgram'], fontsize=14)
plt.show()                              # 上記のグラフをまとめて1つの図に出力

↑

ランダムウォーク（酔歩） †

"np.random..choice([-1,1],M)"を用いて，+1か-1の値がそれぞれ50%の確率で出現する乱数&mathjax{S_i};をM個発生する．

M = 10000                                # 酔歩の回数
np.random.seed(0)                        # 乱数の初期化
S = np.random.choice([-1,1],M)           # +1 or -1 をM個生成し，S(j)に保存 
plt.xlabel('$j$', fontsize=16)
plt.ylabel('$S_j$', fontsize=16)
plt.xlim(0,100)
plt.ylim(-1.5,1.5)
plt.plot(S)
plt.show()

乱数&mathjax{S_i};の分布を確認する．

plt.hist(S, bins=20,normed=True)         # Rに格納された乱数の分布を棒グラフにする．さらに総面積が1になるよう規格化
plt.xlabel('$S$', fontsize=16)
plt.ylabel('$P(S)$', fontsize=16)
plt.legend(['histgram'], fontsize=14)
plt.title(r'Binomial Distribution: p=%f, M=%i' % (p,M))
plt.show()                               # 上記のグラフを出力

乱数&mathjax{S_i};の自己相関関数&mathjax{C(k)=\frac{1}{N}\sum_{i=1}^{N} (S_i - \langle S\rangle)(S_{i+k}-\langle S\rangle)};を計算し，乱数間に相関がないことを確認する．

def auto_correlate(dt):
    cor = np.correlate(dt,dt,mode="full")
    return cor[M-1:]
    return cor
ave = 0.0
corr = np.zeros(M)
corr=auto_correlate(R-ave)/M
plt.plot(corr)
plt.xlim(-100,1000)
plt.ylim(-0.2,1.2)
plt.xlabel('$k$', fontsize=16)
plt.ylabel('$C(k)$', fontsize=16)
plt.show()

原点から出発し，+1か-1のランダムウォーク（酔歩）をM回繰り返した後の位置を&mathjax{R};とする．この試行を独立にL回繰り返し，各回の結果を&mathjax{R_i};とすると，その分布関数は平均&mathjax{\langle R\rangle=M(p-1+q)};，分散&mathjax{\sigma^2=4Mp(1-p)};の正規分布に従うことを確認する．

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の乱数S(j)の合計値のをR(i)に保存
plt.hist(R, bins=20,normed=True)        # Rに格納された乱数の分布を20本の棒グラフにする．さらに総面積が1になるよう規格化
p=0.5
std = np.sqrt(4*M*p*(1-p))              # 正規分布の標準偏差を二項分布に合うように設定
ave = M*(p-1+(1-p))                     # 正規分布の平均値を二項分布に合うように設定
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.xlabel('$R$', fontsize=16)
plt.ylabel('$P((M+R)/2;M)$', fontsize=16)
plt.legend(['Gauss','histgram'], fontsize=14)
plt.title(r'Random walk: p=%f, M=%i' % (p,M))
plt.show()                              # 上記のグラフを出力

↑

宿題 †

ランダムウォーク(+1/-1)をM回繰り返した後の位置&mathjax{R};の分布関数は，本来は二項分布 &mathjax{P(I\equiv(M+R)/2;M)}; で表されるが，&mathjax{M};と&mathjax{p};が大きい時には，平均&mathjax{\langle R\rangle=0};，分散&mathjax{\sigma^2=M};の正規分布で近似することができる（&mathjax{M=10^4};については確認済み）．&mathjax{M=10^2,10^3,10^5};の場合についてもこれを確かめよ．

種々の確率分布関数 †

必要なライブラリの読み込み †

正規分布（Gauss分布） †

二項分布 †

中心極限定理の一例（二項分布→正規分布） †

ランダムウォーク（酔歩） †

宿題 †

Search this site