はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
二項分布の計算とグラフの作成をPythonで行います。
【前の内容】
【他の記事一覧】
【目次】
二項分布
二項分布(Binomial Distribution)の計算と作図を行います。
利用するライブラリを読み込みます。
# 利用するライブラリ import numpy as np from scipy.stats import binom, multinomial # 二項分布, 多項分布 from scipy.special import gamma, loggamma # ガンマ関数, 対数ガンマ関数 import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation
分布の変化をアニメーション(gif画像)で確認するのにmatplotlibライブラリのanimationモジュールを利用します。
定義式の確認
まずは、二項分布の定義式を確認します。
二項分布は、次の式で定義されます。詳しくは「二項分布の定義式 - からっぽのしょこ」を参照してください。
ここで、$x$は成功回数、$\phi$は成功確率($x$が1増える確率)です。
確率変数の値$x$は、$x \in \{0, 1, \cdots, M\}$となります。パラメータ$\phi$は、$\phi \in (0, 1)$を満たす必要があります。また、失敗回数は$M - x$、失敗確率($x$が1増えない確率)は$1 - \phi$で表せます。
この対数をとると、次の式になります。
二項分布の平均と分散は、次の式で計算できます。詳しくは「二項分布の平均と分散の導出：定義式を利用 - からっぽのしょこ」を参照してください。
これらの計算を行いグラフを作成します。
確率の計算
二項分布に従う確率を計算する方法をいくつか確認します。
パラメータを設定します。
# パラメータを指定 phi = 0.3 # 試行回数を指定 M = 10 # 確率変数の値を指定:(x <= M) x = 3
二項分布のパラメータ$0 \leq \phi \leq 1$、試行回数$M$、確率変数がとり得る値$x \in \{0, 1, \cdots, M\}$を指定します。設定した値に従う確率を計算します。
まずは、定義式から確率を計算します。
# 定義式により確率を計算 C = gamma(M + 1) / gamma(M - x + 1) / gamma(x + 1) prob = C * phi**x * (1 - phi)**(M - x) print(prob)
0.2668279319999998
二項分布の定義式
で計算します。$C_{\mathrm{Bin}}$は、二項分布の正規化係数です。
階乗$x!$の計算は、ガンマ関数$\Gamma(x + 1) = x!$に置き換えて計算します。ガンマ関数は、SciPyライブラリのspecialモジュールのgamma()で計算できます。
対数をとった定義式から計算します。
# 対数をとった定義式により確率を計算 log_C = loggamma(M + 1) - loggamma(M - x + 1) - loggamma(x + 1) log_prob = log_C + x * np.log(phi) + (M - x) * np.log(1 - phi) pron = np.exp(log_prob) print(prob, log_prob)
0.2668279319999998 -1.321151277766889
対数をとった定義式
を計算します。
対数をとったガンマ関数はloggamma()で計算できます。
計算結果の指数をとると確率が得られます。
指数と対数の性質より$\exp(\log x) = x$です。
次は、SciPyライブラリのモジュールを使って確率を計算します。
二項分布のモジュールbinomの確率計算メソッドpmf()を使って計算します。
# 二項分布の関数により確率を計算 prob = binom.pmf(k=x, n=M, p=phi) print(prob)
0.26682793200000016
成功回数の引数kにx、試行回数の引数nにM、成功確率の引数pにphiを指定します。
logpmf()だと対数をとった確率を計算します。
# 二項分布の対数をとった関数により確率を計算 log_prob = binom.logpmf(k=x, n=M, p=phi) prob = np.exp(log_prob) print(prob, log_prob)
0.26682793200000016 -1.3211512777668881
計算結果の指数をとると確率が得られます。
以降の計算は、変数とパラメータをベクトル形式にしておく必要があります。
# ベクトルに変換 x_v = np.array([M - x, x]) phi_v = np.array([1.0 - phi, phi])
失敗確率(クラス0の出現確率)$1 - \phi$と成功確率(クラス1の出現確率)$\phi$を持つ配列を作成してphi_vとします。
失敗回数(クラス0の出現頻度)$M - x$と成功回数(クラス1の出現頻度)$x$を持つ配列を作成してx_vとします。
多項分布のモジュールmultinomialの確率計算メソッドpmf()を使って計算します。
# 多項分布の関数により確率を計算 prob = multinomial.pmf(x=x_v, n=M, p=phi_v) print(prob)
0.2668279319999999
出現頻度の引数xにx_v、試行回数の引数nにM、出現確率の引数pにphi_vを指定します。
logpmf()だと対数をとった確率を計算します。
# 多項分布の対数をとった関数により確率を計算 log_prob = multinomial.logpmf(x=x_v, n=M, p=phi_v) prob = np.exp(log_prob) print(prob, log_prob)
0.2668279319999999 -1.321151277766889
計算結果の指数をとると確率が得られます。
統計量の計算
二項分布の平均と分散を計算します。
平均を計算します。
# 平均を計算 E_x = M * phi print(E_x)
3.0
二項分布の平均は、次の式で計算できます。
二項分布のモジュールの平均メソッドmean()でも計算できます。
# 二項分布の関数により平均を計算 print(binom.mean(n=M, p=phi))
3.0
確率計算と同様に引数を指定します。
分散を計算します。
# 分散を計算 V_x = M * phi * (1.0 - phi) print(V_x)
2.0999999999999996
二項分布の分散は、次の式で計算できます。
var()メソッドで分散を計算できます。
# 二項分布の関数により分散を計算 print(binom.var(n=M, p=phi))
2.0999999999999996
こちらも同様に引数を指定します。
グラフの作成
MatplotlibライブラリのPyPlotモジュールを利用して二項分布のグラフを作成します。
二項分布の確率変数$x$がとり得る値ごとの確率を計算します。
# 作図用のxの値を作成 x_vals = np.arange(M + 1) # 分布を計算 probability = binom.pmf(k=x_vals, n=M, p=phi) print(probability[:5])
[0.02824752 0.12106082 0.23347444 0.26682793 0.20012095]
$x$がとり得る値0からMを作成してx_valsとします。
x_valsの各値となる確率を計算します。
二項分布のグラフを作成します。
# 二項分布を作図 plt.figure(figsize=(12, 8)) # 図の設定 plt.bar(x=x_vals, height=probability, color='#00A968') # 棒グラフ #plt.vlines(x=E_x, ymin=0.0, ymax=np.max(probability), color='orange', linestyle='--', label='$E[x]$') # 平均 #plt.vlines(x=E_x - V_x, ymin=0.0, ymax=np.max(probability), color='orange', linestyle=':', label='$E[x] - \sqrt{V[x]}$') # 平均 - 標準偏差 #plt.vlines(x=E_x + V_x, ymin=0.0, ymax=np.max(probability), color='orange', linestyle=':', label='$E[x] + \sqrt{V[x]}$') # 平均 + 標準偏差 plt.xlabel('x') # x軸ラベル plt.ylabel('probability') # y軸ラベル plt.suptitle('Binomial Distribution', fontsize=20) # 図タイトル plt.title('$\phi=' + str(phi) + ', M=' + str(M) + '$', loc='left') # タイトル plt.xticks(ticks=x_vals) # x軸目盛 #plt.legend() # 凡例 plt.grid() # グリッド線 plt.show() # 描画
二項分布のグラフを描画できました。
パラメータとグラフの形状の関係
続いて、パラメータ$\phi$の値を少しずつ変更して、分布の変化をアニメーションで確認します。
・作図コード(クリックで展開)
# 作図用のphiの値を作成 phi_vals = np.arange(start=0.0, stop=1.01, step=0.01) # 図を初期化 fig = plt.figure(figsize=(12, 8)) # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i回目の値を取得 phi = phi_vals[i] # 分布を計算 probability = binom.pmf(k=x_vals, n=M, p=phi) # 二項分布を作図 plt.bar(x=x_vals, height=probability, color='#00A968') # 棒グラフ plt.xlabel('x') # x軸ラベル plt.ylabel('probability') # y軸ラベル plt.suptitle('Binomial Distribution', fontsize=20) # 図タイトル plt.title('$\phi=' + str(np.round(phi, 2)) + ', M=' + str(M) + '$', loc='left') # タイトル plt.xticks(ticks=x_vals) # x軸目盛 plt.grid() # グリッド線 plt.ylim(-0.1, 1.1) # y軸の表示範囲 # gif画像を作成 anime_prob = FuncAnimation(fig, update, frames=len(phi_vals), interval=100) # gif画像を保存 anime_prob.save('Binomial_prob_phi.gif')
$\phi$がとり得る値を作成してphi_valsとします。
phi_valsの値ごとに分布を計算して作図します。
パラメータ$\phi$の値が大きくなるのに従って、成功回数$x$が大きいほど確率が高くなる(山が右に移動する)のを確認できます。
続いて、$\phi$を固定して$M$を変更してみます。
・作図コード(クリックで展開)
# パラメータを指定 phi = 0.3 # 試行回数の最大値を指定 M_max = 100 # 図を初期化 fig = plt.figure(figsize=(12, 8)) # 作図処理を関数として定義 def update(M): # 前フレームのグラフを初期化 plt.cla() # 作図用のxの値を作成 x_vals = np.arange(M + 1) # 分布を計算 probability = binom.pmf(k=x_vals, n=M, p=phi) # 二項分布を作図 plt.bar(x=x_vals, height=probability, color='#00A968') # 棒グラフ plt.xlabel('x') # x軸ラベル plt.ylabel('probability') # y軸ラベル plt.suptitle('Binomial Distribution', fontsize=20) # 図タイトル plt.title('$\phi=' + str(np.round(phi, 2)) + ', M=' + str(M) + '$', loc='left') # タイトル #plt.xticks(ticks=x_vals) # x軸目盛 plt.grid() # グリッド線 plt.ylim(-0.1, 1.1) # y軸の表示範囲 # gif画像を作成 anime_prob = FuncAnimation(fig, update, frames=M_max, interval=100) # gif画像を保存 anime_prob.save('Binomial_prob_M.gif')
試行回数$M$の最大値(フレーム数)をM_maxに指定します。
0からM_max - 1まで順番にMの値が更新され、x_valsの要素も更新されます。
更新されたx_valsに対して分布を計算して作図します。
$M$が増えるのに従って、成功回数$x$が大きいほど確率が高くなる(山が右に移動する)のを確認できます。ただし、$x$がとり得る値の範囲x_valsも広がっていくため、x_vals全体における山の位置は変わりません。
乱数の生成
二項分布の乱数を生成してヒストグラムを確認します。
パラメータを指定して、二項分布に従う乱数を生成します。
# パラメータを指定 phi = 0.3 # 試行回数を指定 M = 10 # データ数を指定 N = 1000 # 二項分布に従う乱数を生成 x_n = np.random.binomial(n=M, p=phi, size=N) # 確認 print(x_n[:5])
[1 6 2 1 1]
二項分布の乱数生成関数np.random.binomial()の試行回数の引数nにM、成功確率の引数pにphi、データ数(サンプルサイズ)の引数sizeにNを指定します。
サンプルの値を集計します。
# 乱数を集計 frequency = np.array([np.sum(x_n == m) for m in range(M + 1)]) print(frequency[:5])
[ 30 113 248 277 179]
x_nに含まれる値がmの要素数は、np.sum(x_n == m)で得られます。
リスト内包表記で、0からMまでカウントして配列に格納します。
ヒストグラムを作成します。
# サンプルのヒストグラムを作成 plt.figure(figsize=(12, 8)) # 図の設定 plt.bar(x=x_vals, height=frequency, color='#00A968') # ヒストグラム plt.xlabel('x') # x軸ラベル plt.ylabel('frequency') # y軸ラベル plt.suptitle('Binomial Distribution', fontsize=20) # 図タイトル plt.title('$\phi=' + str(phi) + ', N=' + str(N) + '=(' + ', '.join([str(f) for f in frequency]) + ')$', loc='left') # タイトル plt.xticks(ticks=x_vals) # x軸目盛 plt.grid() # グリッド線 plt.show() # 描画
構成比を分布と重ねて描画します。
# サンプルの構成比を作図 plt.figure(figsize=(12, 8)) # 図の設定 plt.bar(x=x_vals, height=probability, color='white', edgecolor='green', linestyle='--') # 分布 plt.bar(x=x_vals, height=frequency / N, color='#00A968', alpha=0.8) # 構成比 plt.xlabel('x') # x軸ラベル plt.ylabel('proportion') # y軸ラベル plt.suptitle('Binomial Distribution', fontsize=20) # 図タイトル plt.title('$\phi=' + str(phi) + ', N=' + str(N) + '=(' + ', '.join([str(f) for f in frequency]) + ')$', loc='left') # タイトル plt.xticks(ticks=x_vals) # x軸目盛 plt.grid() # グリッド線 plt.show() # 描画
頻度frequencyをデータ数Nで割り、各値の構成比を計算します。
データ数が十分に増えると分布のグラフに形が近づきます。
サンプルサイズとヒストグラムの変化をアニメーションで確認します。乱数を1つずつ取り出して作図します。
ヒストグラムのアニメーションを作成します。
・作図コード(クリックで展開)
# フレーム数を指定 N = 100 # 図を初期化 fig = plt.figure(figsize=(12, 8)) # 頻度の最大値を取得 y_max = np.max([np.sum(x_n[:N] == m) for m in range(M + 1)]) # 作図処理を関数として定義 def update(n): # 前フレームのグラフを初期化 plt.cla() # n個の乱数を集計 frequency = np.array([np.sum(x_n[:(n+1)] == m) for m in range(M + 1)]) # サンプルのヒストグラムを作成 plt.bar(x=x_vals, height=frequency, color='#00A968', zorder=1) # ヒストグラム plt.scatter(x=x_n[n], y=0.0, color='orange', s=100, zorder=2) # サンプル plt.xlabel('x') # x軸ラベル plt.ylabel('frequency') # y軸ラベル plt.suptitle('Binomial Distribution', fontsize=20) # 図タイトル plt.title('$\phi=' + str(phi) + ', N=' + str(n) + '=(' + ', '.join([str(f) for f in frequency]) + ')$', loc='left') # タイトル plt.xticks(ticks=x_vals) # x軸目盛 plt.grid() # グリッド線 plt.ylim(-1.0, y_max + 1.0) # y軸の表示範囲 # gif画像を作成 anime_freq = FuncAnimation(fig, update, frames=N, interval=100) # gif画像を保存 anime_freq.save('Binomial_freq.gif')
構成比のアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(9, 6)) # 作図処理を関数として定義 def update(n): # 前フレームのグラフを初期化 plt.cla() # n個の乱数を集計 frequency = np.array([np.sum(x_n[:(n+1)] == m) for m in range(M + 1)]) # サンプルの構成比を作成 plt.bar(x=x_vals, height=probability, color='white', edgecolor='green', linestyle='--', zorder=1) # 分布 plt.bar(x=x_vals, height=frequency / (n + 1), color='#00A968', alpha=0.8, zorder=2) # 構成比 plt.scatter(x=x_n[n], y=0.0, color='orange', s=100, zorder=3) # サンプル plt.xlabel('x') # x軸ラベル plt.ylabel('proportion') # y軸ラベル plt.suptitle('Binomial Distribution', fontsize=20) # 図タイトル plt.title('$\phi=' + str(phi) + ', N=' + str(n) + '=(' + ', '.join([str(f) for f in frequency]) + ')$', loc='left') # タイトル plt.xticks(ticks=x_vals) # x軸目盛 plt.grid() # グリッド線 plt.ylim(-0.1, 1.1) # y軸の表示範囲 # gif画像を作成 anime_prop = FuncAnimation(fig, update, frames=N, interval=100) # gif画像を保存 anime_prop.save('Binomial_prop.gif')
サンプルが増えるに従って、真の分布に近付いていくのを確認できます。
以上で、二項分布を確認できました。
参考文献
- 岩田具治『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社,2015年.
おわりに
投稿日に放送開始されたこのドラマを観てね。
ハロー！ワールド
【次の内容】
