はじめに
機械学習で登場する確率分布について色々な角度から理解したいシリーズです。
カテゴリ分布の計算とグラフの作成をPythonで行います。
【前の内容】
【他の記事一覧】
【目次】
カテゴリ分布
カテゴリ分布(Categorical Distribution)の計算と作図を行います。
利用するライブラリを読み込みます。
# 利用するライブラリ import numpy as np from scipy.stats import multinomial # 多項分布 import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation
分布の変化をアニメーション(gif画像)で確認するのにmatplotlibライブラリのanimationモジュールを利用します。
定義式の確認
まずは、カテゴリ分布の定義式を確認します。
カテゴリ分布は、次の式で定義されます。詳しくは「カテゴリ分布の定義式 - からっぽのしょこ」を参照してください。
ここで、$x_v$はクラス$v$が出現した回数、$\phi_v$はクラス$v$の出現確率です。
確率変数の値$\mathbf{x} = (x_1, \cdots, x_V)$はone-hotベクトルで、$x_v \in \{0, 1\}$、$\sum_{v=1}^V x_v = 1$となります。パラメータ$\boldsymbol{\phi} = (\phi_1, \cdots, \phi_V)$は、$\phi_v \in (0, 1)$、$\sum_{v=1}^V \phi_v = 1$を満たす必要があります。
この式の対数をとると、次の式になります。
カテゴリ分布のクラス$v$における平均と分散は、次の式で計算できます。詳しくは「カテゴリ分布の平均と分散の導出 - からっぽのしょこ」を参照してください。
これらの計算を行いグラフを作成します。
統計量の計算
カテゴリ分布の平均と分散を計算します。
クラス$v$の平均を計算します。
# クラスvの平均を計算 E_x = phi_v[v] print(E_x)
0.4
カテゴリ分布の平均は、次の式で計算できます。
クラス$v$の分散を計算します。
# クラスvの分散を計算 V_x = phi_v[v] * (1.0 - phi_v[v]) print(V_x)
0.24
カテゴリ分布の分散は、次の式で計算できます。
グラフの作成
MatplotlibライブラリのPyPlotモジュールを利用してカテゴリ分布のグラフを作成します。
カテゴリ分布の確率変数がとり得るクラス$v$と対応する確率を作成します。
# パラメータを指定 phi_v = np.array([0.2, 0.4, 0.1, 0.3]) # クラス数を取得 V = len(phi_v) # 作図用のクラス番号を作成 v_vals = np.arange(1, V + 1) # 分布を計算 probability = phi_v.copy()
$\mathbf{x}$によって表されるクラス番号1からVをv_vals、各クラスに対応する確率phi_vをprobabilityとします。
multinomial.pmf()でも計算できます。
# 多項分布の関数により分布を計算 probability = multinomial.pmf(x=np.identity(V), n=1, p=phi_v) print(probability)
[0.2 0.4 0.1 0.3]
単位行列np.identity(V)を使うことで、各クラスの確率phi_vを取り出せます。
カテゴリ分布のグラフを作成します。
# カテゴリ分布を作図 plt.figure(figsize=(12, 8)) # 図の設定 plt.bar(x=v_vals, height=probability, color='#00A968') # 棒グラフ plt.xlabel('v') # x軸ラベル plt.ylabel('probability') # y軸ラベル plt.suptitle('Categorical Distribution', fontsize=20) # 図タイトル plt.title('$\phi=(' + ', '.join([str(phi) for phi in phi_v]) + ')$', loc='left') # タイトル plt.xticks(ticks=v_vals) # x軸目盛 plt.grid() # グリッド線 plt.show() # 描画
パラメータの値そのままですが、これがカテゴリ分布のグラフです。
パラメータと分布の形状の関係
続いて、パラメータ$\phi_v$の値を少しずつ変更して、分布の変化をアニメーションで確認します。
・作図コード(クリックで展開)
# 作図用のphi_1の値を作成 phi_vals = np.arange(start=0.0, stop=1.0, step=0.01) # クラス数を指定 V = 3 # 作図用のクラス番号を作成 v_vals = np.arange(1, V + 1) # 図を初期化 fig = plt.figure(figsize=(12, 8)) # 作図処理を関数として定義 def update(i): # 前フレームのグラフを初期化 plt.cla() # i回目のphi_1の値を取得 phi_1 = phi_vals[i] # phi_1以外の割り当てを指定 phi_v = np.array([phi_1, (1.0 - phi_1) * 0.6, (1.0 - phi_1) * 0.4]) # カテゴリ分布を作図 plt.bar(x=v_vals, height=phi_v, color='#00A968', zorder=1) # 棒グラフ plt.xlabel('v') # x軸ラベル plt.ylabel('Probability') # y軸ラベル plt.suptitle('Categorical Distribution', fontsize=20) # 図タイトル plt.title('$\phi=(' + ', '.join([str(phi) for phi in np.round(phi_v, 2)]) + ')$', loc='left') # タイトル plt.xticks(ticks=v_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('Categorical_prob.gif')
$\phi_1$がとり得る値を作成してphi_valsとします。
phi_valsの値ごとにphi_vを更新して作図します。$\phi_1$以外の確率の和$1 - \phi_1$を他のクラスの確率として割り振ります。
(定義のままですが、)パラメータ$\phi_1$の値が大きいほど、$x_1 = 1$となる確率が高く、他のクラスとなる確率が低いを確認できます。
乱数の生成
カテゴリ分布の乱数を生成してヒストグラムを確認します。
パラメータを指定して、カテゴリ分布に従う乱数を生成します。
# パラメータを指定 phi_v = np.array([0.2, 0.4, 0.1, 0.3]) # クラス数を取得 V = len(phi_v) # 作図用のクラス番号を作成 v_vals = np.arange(1, V + 1) # データ数を指定 N = 1000 # カテゴリ分布に従う乱数を生成 x_nv = np.random.multinomial(n=1, pvals=phi_v, size=N) print(x_nv[:5])
[[0 1 0 0]
[1 0 0 0]
[0 1 0 0]
[0 0 0 1]
[0 1 0 0]]
多項分布の乱数生成関数np.random.multinomial()の試行回数の引数nに1を指定することで、カテゴリ分布に従う乱数を生成できます。
成功確率の引数pvalsにphi、データ数(サンプルサイズ)の引数sizeにNを指定します。
各データに割り当てられたクラス番号を抽出します。
# クラス番号を抽出 x_n = np.where(x_nv == 1)[1] print(x_n[:5])
[1 0 1 3 1]
np.where()を使って、x_nvの各行から値が1の列番号を抽出できます。ここではx_nは使いません。
サンプルの値を集計します。
# 乱数を集計 frequency = np.sum(x_nv, axis=0) print(frequency)
[198 424 80 298]
x_nvの各列に含まれる1の要素数は、np.sum(x_n == 0, axis=0)で得られます。
ヒストグラムを作成します。
# サンプルのヒストグラムを作成 plt.figure(figsize=(12, 8)) # 図の設定 plt.bar(x=v_vals, height=frequency, color='#00A968') # ヒストグラム plt.xlabel('v') # x軸ラベル plt.ylabel('frequency') # y軸ラベル plt.suptitle('Bernoulli Distribution', fontsize=20) # 図タイトル plt.title('$\phi=(' + ', '.join([str(phi) for phi in phi_v]) + ')' + ', N=' + str(N) +'=(' + ', '.join([str(f) for f in frequency]) + ')$', loc='left') # タイトル plt.xticks(ticks=v_vals) # x軸目盛 plt.grid() # グリッド線 plt.show() # 描画
構成比を分布と重ねて描画します。
# サンプルの構成比を作図 plt.figure(figsize=(12, 8)) # 図の設定 plt.bar(x=v_vals, height=probability, color='white', edgecolor='green', linestyle='--') # 分布 plt.bar(x=v_vals, height=frequency / N, color='#00A968', alpha=0.8) # 構成比 plt.xlabel('v') # x軸ラベル plt.ylabel('proportion') # y軸ラベル plt.suptitle('Bernoulli Distribution', fontsize=20) # 図タイトル plt.title('$\phi=(' + ', '.join([str(phi) for phi in phi_v]) + ')' + ', N=' + str(N) +'=(' + ', '.join([str(f) for f in frequency]) + ')$', loc='left') # タイトル plt.xticks(ticks=v_vals) # x軸目盛 plt.grid() # グリッド線 plt.show() # 描画
各クラスの頻度frequencyをデータ数Nで割り、0と1の構成比を計算します。
データ数が十分に増えると分布に形が近づきます。
サンプルサイズとヒストグラムの変化をアニメーションで確認します。乱数を1つずつ取り出して作図します。
・作図コード(クリックで展開)
ヒストグラムのアニメーションを作成します。
# フレーム数を指定 N = 100 # 図を初期化 fig = plt.figure(figsize=(12, 8)) # 頻度の最大値を取得 y_max = np.max(np.sum(x_nv[:N], axis=0)) # 作図処理を関数として定義 def update(n): # 前フレームのグラフを初期化 plt.cla() # n個の乱数を集計 frequency = np.sum(x_nv[:(n+1)], axis=0) # n番目の乱数のクラスを取得 x_val = np.where(x_nv[n] == 1)[0][0] + 1 # サンプルのヒストグラムを作成 plt.bar(x=v_vals, height=frequency, color='#00A968', zorder=1) # ヒストグラム plt.scatter(x=x_val, y=0.0, color='orange', s=100, zorder=2) # サンプル plt.xlabel('v') # x軸ラベル plt.ylabel('frequency') # y軸ラベル plt.suptitle('Categorical Distribution', fontsize=20) # 図タイトル plt.title('$\phi=(' + ', '.join([str(phi) for phi in phi_v]) + ')' + ', N=' + str(n + 1) +'=(' + ', '.join([str(f) for f in frequency]) + ')$', loc='left') # タイトル plt.xticks(ticks=v_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('Categorical_freq.gif')
構成比のアニメーションを作成します。
# 図を初期化 fig = plt.figure(figsize=(12, 8)) # 作図処理を関数として定義 def update(n): # 前フレームのグラフを初期化 plt.cla() # n個の乱数を集計 frequency = np.sum(x_nv[:(n+1)], axis=0) # n番目の乱数のクラスを取得 x_val = np.where(x_nv[n] == 1)[0][0] + 1 # サンプルの構成比を作成 plt.bar(x=v_vals, height=probability, color='white', edgecolor='green', linestyle='--', zorder=1) # 分布 plt.bar(x=v_vals, height=frequency / (n + 1), color='#00A968', alpha=0.8, zorder=2) # 構成比 plt.scatter(x=x_val, y=0.0, color='orange', s=100, zorder=3) # サンプル plt.xlabel('v') # x軸ラベル plt.ylabel('proportion') # y軸ラベル plt.suptitle('Categorical Distribution', fontsize=20) # 図タイトル plt.title('$\phi=(' + ', '.join([str(phi) for phi in phi_v]) + ')' + ', N=' + str(n + 1) +'=(' + ', '.join([str(f) for f in frequency]) + ')$', loc='left') # タイトル plt.xticks(ticks=v_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('Categorical_prop.gif')
サンプルが増えるに従って、真の分布に近付いていくのを確認できます。
以上で、カテゴリ分布を確認できました。
参考文献
- 岩田具治『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社,2015年.
- 須山敦志『ベイズ推論による機械学習入門』(機械学習スタートアップシリーズ)杉山将監修,講談社,2017年.
おわりに
ここまでは準備運動みたいなもので、次から変則的な作図になります。
【次の内容】
