はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
ベルヌーイ分布の計算とグラフの作成をPythonで行います。
【他の記事一覧】
【この記事の内容】
ベルヌーイ分布
ベルヌーイ分布(Bernoulli Distribution)の計算と作図を行います。
利用するライブラリを読み込みます。
# 利用するライブラリ import numpy as np from scipy.stats import bernoulli, binom, multinomial # ベルヌーイ分布, 二項分布, 多項分布 import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation
分布の変化をアニメーション(gif画像)で確認するのにmatplotlib
ライブラリのanimation
モジュールを利用します。
定義式の確認
まずは、ベルヌーイ分布の定義式を確認します。
ベルヌーイ分布は、次の式で定義されます。詳しくは「ベルヌーイ分布の定義式 - からっぽのしょこ」を参照してください。
ここで、$x$は成功・失敗を表す値、$\phi$は成功確率($x = 1$となる確率)です。
確率変数の値$x$は、$x \in \{0, 1\}$となり、$x = 1$が成功・$x = 0$が失敗を表します。パラメータ$\phi$は、$\phi \in (0, 1)$を満たす必要があります。また、失敗確率($x = 0$となる確率)は$1 - \phi$で表せます。
この式の対数をとると、次の式になります。
ベルヌーイ分布の平均と分散は、次の式で計算できます。詳しくは「ベルヌーイ分布の平均と分散の導出 - からっぽのしょこ」を参照してください。
これらの計算を行いグラフを作成します。
確率の計算
ベルヌーイ分布に従う確率を計算する方法をいくつか確認します。
パラメータを設定します。
# パラメータを指定 phi = 0.3 # 確率変数の値を指定 x = 1.0
ベルヌーイ分布のパラメータ$0 \leq \phi \leq 1$と確率変数がとり得る値$x \in {0, 1}$を指定します。設定した値に従う確率を計算します。
まずは、定義式から確率を計算します。
# 定義式により確率を計算 prob = phi**x * (1.0 - phi)**(1.0 - x) print(prob)
0.3
ベルヌーイ分布の定義式
で計算します。
対数をとった定義式から計算します。
# 対数をとった定義式により確率を計算 log_prob = x * np.log(phi) + (1.0 - x) * np.log(1.0 - phi) prob = np.exp(log_prob) print(prob, log_prob)
0.3 -1.2039728043259361
対数をとった定義式
を計算します。計算結果の指数をとると確率が得られます。
指数と対数の性質より$\exp(\log x) = x$です。
次は、SciPy
ライブラリのモジュールを使って確率を計算します。
ベルヌーイ分布のモジュールbernoulli
の確率計算メソッドpmf()
を使って計算します。
# ベルヌーイ分布の関数により確率を計算 prob = bernoulli.pmf(k=x, p=phi) print(prob)
0.3
引数k
にx
、成功確率の引数p
にphi
を指定します。
logpmf()
だと対数をとった確率を計算します。
# ベルヌーイ分布の対数をとった関数により確率を計算 log_prob = bernoulli.logpmf(k=x, p=phi) prob = np.exp(log_prob) print(prob, log_prob)
0.3 -1.2039728043259361
計算結果の指数をとると確率が得られます。
二項分布のモジュールbinom
の確率計算メソッドpmf()
を使って計算します。
# 二項分布の関数により確率を計算 prob = binom.pmf(k=x, n=1, p=phi) print(prob)
0.3
試行回数の引数n
に1
を指定することで、ベルヌーイ分布の確率を計算できます。
成功回数の引数k
にx
、成功確率の引数p
にphi
を指定します。
logpmf()
だと対数をとった確率を計算します。
# 二項分布の対数をとった関数により確率を計算 log_prob = binom.logpmf(k=x, n=1, p=phi) prob = np.exp(log_prob) print(prob, log_prob)
0.3 -1.2039728043259361
計算結果の指数をとると確率が得られます。
以降の計算は、変数とパラメータをベクトルに変換しておく必要があります。
# ベクトルに変換 x_v = np.array([1 - x, x]) phi_v = np.array([1.0 - phi, phi])
$x = 0$となる確率$1 - \phi$と、$x = 1$となる確率$\phi$を持つ配列を作成してphi_v
とします。
$x = 0$のとき$1 - x = 1$、$x = 1$のとき$1 - x = 0$となるのを利用して、one-hotベクトル(1-of-K符号化法)の$x$を作成してx_v
とします。
多項分布のモジュールmultinomial
の確率計算メソッドpmf()
を使って計算します。
# 多項分布の関数により確率を計算 prob = multinomial.pmf(x=x_v, n=1, p=phi_v) print(prob)
0.30000000000000004
二項分布と同様に試行回数の引数をn=1
として、出現頻度の引数x
にx
、出現確率の引数p
にphi_v
を指定します。
logpmf()
だと対数をとった確率を計算します。
# 多項分布の対数をとった関数により確率を計算 log_prob = multinomial.logpmf(x=x_v, n=1, p=phi_v) prob = np.exp(log_prob) print(prob, log_prob)
0.30000000000000004 -1.203972804325936
計算結果の指数をとると確率が得られます。
最後に、スライス機能を使って確率を取り出します。
# インデックスにより確率を抽出 prob = phi_v[np.int(x)] print(prob)
0.3
x
が0
となる確率はphi_v
の0番目の要素で、1
となる確率はphi_v
の1番目の要素なので、x
の値に対応するパラメータのインデックスはx
になります。
統計量の計算
ベルヌーイ分布の平均と分散を計算します。
平均を計算します。
# 定義式により平均を計算 E_x = phi print(E_x)
0.3
ベルヌーイ分布の平均は、次の式で計算できます。
各分布のモジュールの平均メソッドmean()
でも計算できます。
# ベルヌーイ分布の関数により平均を計算 print(bernoulli.mean(p=phi)) # 二項分布の関数により平均を計算 print(binom.mean(n=1, p=phi))
0.3
0.3
確率計算と同様に引数を指定します。
分散を計算します。
# 分散を計算 V_x = phi * (1.0 - phi) print(V_x)
0.21
ベルヌーイ分布の分散は、次の式で計算できます。
var()
メソッドで分散を計算できます。
# ベルヌーイ分布の関数により分散を計算 print(bernoulli.var(p=phi)) # 二項分布の関数により分散を計算 print(binom.var(n=1, p=phi))
0.21
0.21
こちらも同様に引数を指定します。
グラフの作成
Matplotlib
ライブラリのPyPlot
モジュールを利用してベルヌーイ分布のグラフを作成します。
作図用の配列を作成します。
# パラメータを指定 phi = 0.3 # 作図用の値を作成 x_vals = np.array([0.0, 1.0]) probability = np.array([1.0 - phi, phi]) print(x_vals) print(probability)
[0. 1.]
[0.7 0.3]
$x$がとり得る値0, 1
と、それに対応する確率1 - phi, phi
を配列に格納します。
ベルヌーイ分布のグラフを作成します。
# ベルヌーイ分布を作図 plt.figure(figsize=(9, 8)) # 図の設定 plt.bar(x_vals, 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('Bernoulli Distribution', fontsize=20) # 図タイトル plt.title('$\phi=' + str(phi) + '$', loc='left') # タイトル plt.xticks(ticks=[0, 1]) # 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=(9, 8)) # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i回目の値を取得 phi = phi_vals[i] # ベルヌーイ分布を作図 plt.bar([0.0, 1.0], [1.0 - phi, phi], color='#00A968') # 棒グラフ plt.xlabel('x') # x軸ラベル plt.ylabel('probability') # y軸ラベル plt.suptitle('Bernoulli Distribution', fontsize=20) # 図タイトル plt.title('$\phi=' + str(np.round(phi, 2)) + '$', loc='left') # タイトル plt.xticks(ticks=[0, 1]) # 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('Bernoulli_prob.gif')
$\phi$がとり得る値を作成してphi_vals
とします。
phi_vals
の値ごとに作図します。
(定義のままですが、)パラメータ$\phi$の値が大きくなるほど、$x = 0$となる確率が下がり、$x = 1$となる確率が上がるのを確認できます。
乱数の生成
ベルヌーイ分布の乱数を生成してヒストグラムを確認します。
パラメータを指定して、ベルヌーイ分布に従う乱数を生成します。
# パラメータを指定 phi = 0.3 # データ数を指定 N = 1000 # ベルヌーイ分布に従う乱数を生成 x_n = np.random.binomial(n=1, p=phi, size=N) print(x_n[:5])
[0 0 0 0 0]
二項分布の乱数生成関数np.random.binomial()
の試行回数の引数n
に1
を指定することで、ベルヌーイ分布に従う乱数を生成できます。
データ数(サンプルサイズ)の引数size
にN
、成功確率の引数p
にphi
を指定します。
サンプルの値を集計します。
# 乱数を集計 frequency = np.array([np.sum(x_n == 0), np.sum(x_n == 1)]) print(frequency)
[685 315]
x_n
に含まれる0
の要素数は、np.sum(x_n == 0)
で得られます。値が1
の要素数も同様に求めて、配列に格納します。
また、得られた頻度frequency
をデータ数N
で割り、0
と1
の構成比を計算します。
ヒストグラムを作成します。
# サンプルのヒストグラムを作成 plt.figure(figsize=(9, 8)) # 図の設定 plt.bar(x_vals, frequency, color='#00A968') # ヒストグラム plt.xlabel('x') # x軸ラベル plt.ylabel('frequency') # y軸ラベル plt.suptitle('Bernoulli Distribution', fontsize=20) # 図タイトル plt.title('$\phi=' + str(phi) + ', N=' + str(N) + '=(' + str(np.sum(x_n == 0)) + ', ' + str(np.sum(x_n == 1)) + ')$', loc='left') # タイトル plt.xticks(ticks=[0, 1]) # x軸目盛 plt.grid() # グリッド線 plt.show() # 描画
構成比を分布と重ねて描画します。
# サンプルの構成比を作図 plt.figure(figsize=(9, 8)) # 図の設定 plt.bar(x_vals, probability, color='white', edgecolor='green', linestyle='--') # 分布 plt.bar(x_vals, frequency / N, color='#00A968', alpha=0.8) # 構成比 plt.xlabel('x') # x軸ラベル plt.ylabel('proportion') # y軸ラベル plt.suptitle('Bernoulli Distribution', fontsize=20) # 図タイトル plt.title('$\phi=' + str(phi) + ', N=' + str(N) + '=(' + str(np.sum(x_n == 0)) + ', ' + str(np.sum(x_n == 1)) + ')$', loc='left') # タイトル plt.xticks(ticks=[0, 1]) # x軸目盛 plt.grid() # グリッド線 plt.show() # 描画
データ数が十分に増えると分布に形が近づきます。
サンプルサイズとヒストグラムの変化をアニメーションで確認します。
・作図コード(クリックで展開)
ヒストグラムのアニメーションを作成します。
乱数を1つずつ取り出して作図します。
# フレーム数を指定 N = 100 # 図を初期化 fig = plt.figure(figsize=(9, 8)) # 頻度の最大値を取得 y_max = np.max([np.sum(x_n[:N] == 0), np.sum(x_n[:N] == 1)]) # 作図処理を関数として定義 def update(n): # 前フレームのグラフを初期化 plt.cla() # n個の乱数を集計 frequency = np.array([np.sum(x_n[:(n+1)] == 0), np.sum(x_n[:(n+1)] == 1)]) # サンプルのヒストグラムを作成 plt.bar(x_vals, frequency, color='#00A968', zorder=1) # ヒストグラム plt.scatter(x_n[n], 0.0, color='orange', s=100, zorder=2) # サンプル plt.xlabel('x') # x軸ラベル plt.ylabel('frequency') # y軸ラベル plt.suptitle('Bernoulli Distribution', fontsize=20) # 図タイトル plt.title('$\phi=' + str(phi) + ', N=' + str(n) + '=(' + str(frequency[0]) + ', ' + str(frequency[1]) + ')$', loc='left') # タイトル plt.xticks(ticks=[0, 1]) # x軸目盛 plt.grid() # グリッド線 plt.ylim(-1.0, y_max + 1.0) # y軸の表示範囲 # gif画像を作成 anime_hist = FuncAnimation(fig, update, frames=N, interval=100) # gif画像を保存 anime_hist.save('Bernoulli_hist.gif')
構成比のアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(9, 8)) # 作図処理を関数として定義 def update(n): # 前フレームのグラフを初期化 plt.cla() # n個の乱数を集計 frequency = np.array([np.sum(x_n[:(n+1)] == 0), np.sum(x_n[:(n+1)] == 1)]) # サンプルの構成比を作成 plt.bar(x_vals, probability, color='white', edgecolor='green', linestyle='--', zorder=1) # 分布 plt.bar(x_vals, frequency / (n + 1), color='#00A968', alpha=0.8, zorder=2) # 構成比 plt.scatter(x_n[n], 0.0, color='orange', s=100, zorder=3) # サンプル plt.xlabel('x') # x軸ラベル plt.ylabel('proportion') # y軸ラベル plt.suptitle('Bernoulli Distribution', fontsize=20) # 図タイトル plt.title('$\phi=' + str(phi) + ', N=' + str(n) + '=(' + str(frequency[0]) + ', ' + str(frequency[1]) + ')$', loc='left') # タイトル plt.xticks(ticks=[0, 1]) # x軸目盛 plt.grid() # グリッド線 plt.ylim(-0.1, 1.1) # y軸の表示範囲 # gif画像を作成 anime_prop = FuncAnimation(fig, update, frames=N, interval=100) # gif画像を保存 anime_prop.save('Bernoulli_prop.gif')
サンプルが増えるに従って、真の分布に近付いていくのを確認できます。
以上で、ベルヌーイ分布を確認できました。
参考文献
- 岩田具治『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社,2015年.
おわりに
PyPlot
とggplot2
とでlinestyle
とlinetype
みたいな引数名の微妙な違いはしょうがないとして、SciPy
とstats
とでn
とsize
が逆なのは厄介。
【次の内容】