はじめに

　機械学習で登場する確率分布について色々な角度から理解したいシリーズです。

　ベルヌーイ分布の計算とグラフの作成をPythonで行います。

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

はじめに
ベルヌーイ分布
参考文献
おわりに

ベルヌーイ分布

　ベルヌーイ分布(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モジュールを利用します。

定義式の確認

　まずは、ベルヌーイ分布の定義式を確認します。

　ベルヌーイ分布は、次の式で定義されます。詳しくは「ベルヌーイ分布の定義式 - からっぽのしょこ」を参照してください。

$$ \mathrm{Bern}(x | \phi) = \phi^x (1 - \phi)^{1 - x} $$

　ここで、$x$は成功・失敗を表す値、$\phi$は成功確率($x = 1$となる確率)です。
　確率変数の値$x$は、$x \in \{0, 1\}$となり、$x = 1$が成功・$x = 0$が失敗を表します。パラメータ$\phi$は、$\phi \in (0, 1)$を満たす必要があります。また、失敗確率($x = 0$となる確率)は$1 - \phi$で表せます。

　この式の対数をとると、次の式になります。

$$ \log \mathrm{Bern}(x | \phi) = x \log \phi + (1 - x) \log (1 - \phi) $$

　ベルヌーイ分布の平均と分散は、次の式で計算できます。詳しくは「ベルヌーイ分布の平均と分散の導出 - からっぽのしょこ」を参照してください。

$$ \begin{aligned} \mathbb{E}[x] &= \phi \\ \mathbb{V}[x] &= \phi (1 - \phi) \end{aligned} $$

　これらの計算を行いグラフを作成します。

確率の計算

　ベルヌーイ分布に従う確率を計算する方法をいくつか確認します。

　パラメータを設定します。

# パラメータを指定
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

　ベルヌーイ分布の定義式

$$ p(x | \phi) = \phi^x (1 - \phi)^{1 - x} $$

で計算します。
　対数をとった定義式から計算します。

# 対数をとった定義式により確率を計算
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

　対数をとった定義式

$$ \log p(x | \phi) = x \log \phi + (1 - x) \log (1 - \phi) $$

を計算します。計算結果の指数をとると確率が得られます。

$$ p(x | \phi) = \exp \Bigr( \log p(x | \phi) \Bigr) $$

　指数と対数の性質より$\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

　ベルヌーイ分布の平均は、次の式で計算できます。

$$ \mathbb{E}[x] = \phi $$

　各分布のモジュールの平均メソッド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

　ベルヌーイ分布の分散は、次の式で計算できます。

$$ \mathbb{V}[x] = \phi (1 - \phi) $$

　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が逆なのは厄介。

【次の内容】

からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

【Python】ベルヌーイ分布の作図

はじめに

ベルヌーイ分布

定義式の確認

確率の計算

統計量の計算

グラフの作成

パラメータとグラフの形状の関係

乱数の生成

参考文献

おわりに