【Python】ディガンマ関数の作図

はじめに

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

　ディガンマ関数のグラフをPythonで作成します。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

ディガンマ関数の作図

　ディガンマ関数(Digamma Function)の計算を確認して、グラフを作成します。

　利用するライブラリを読み込みます。

# 利用するライブラリ
import numpy as np
from scipy.special import digamma # ディガンマ関数
import matplotlib.pyplot as plt

ディガンマ関数の定義

　まずは、ディガンマ関数の定義を確認します。

　ディガンマ関数は、ガンマ関数の対数の微分で定義されます。ガンマ関数については「ガンマ関数の性質の導出 - からっぽのしょこ」を参照してください。

$$ \psi(x) = \frac{d}{d x} \log \Gamma(x) = \frac{1}{\Gamma(x)} \frac{d \Gamma(x)}{d x} $$

　ここで、$\Gamma(x)$はガンマ関数です。対数関数の微分は$\frac{d \log x}{d x} = \frac{1}{x}$です。

ガンマ関数の計算

　まずは、ディガンマ関数の計算方法を確認します。

　SciPyライブラリのspecialモジュールのディガンマ関数digamma()を使って計算します。

# 変数の値を指定
x = 4.0

# ディガンマ関数の計算
y = digamma(x)
print(y)

1.2561176684318003

　$x$が0または負の整数の場合は計算できません。

# ディガンマ関数を計算:(計算できない)
print(digamma(0.0))
print(digamma(-3.0))

-inf
nan

　ディガンマ関数の計算を確認できました。次は、ディガンマ関数を可視化します。

グラフの作成

　MatplotlibライブラリのPyPlotモジュールを利用して、ディガンマ関数のグラフを作成します。

非負の値の場合

　まずは、0より大きい実数$x$を範囲として、ディガンマ関数を計算します。

# x軸の値を指定:(x > 0)
x_vals = np.arange(start=0.01, stop = 5.0, step=0.01)
print(x_vals[:5])

# ディガンマ関数の計算
y_vals = digamma(x_vals)
print(y_vals[:5])

[0.01 0.02 0.03 0.04 0.05]
[-100.56088546  -50.54478931  -33.86225442  -25.51327488  -20.49784499]

　x軸の値を作成してx_valsとします。
　作成したx_valsを用いて、digamma()でディガンマ関数の計算をします。

　ディガンマ関数のグラフを作成します。

# ディガンマ関数を作図
plt.figure(figsize=(12, 9)) # 図の設定
plt.plot(x_vals, y_vals, color='orange') # 折れ線グラフ
plt.xlabel('x') # x軸ラベル
plt.ylabel('$\psi(x)$') # y軸ラベル
plt.suptitle('Digamma Function', fontsize=20) # 全体のタイトル
plt.grid() # グリッド線
plt.show() # 描画

ディガンマ関数のグラフ

　$x = 0$の垂線に漸近しているのが分かります。

負の値を含む場合

　続いて、負の実数も含めて、ディガンマ関数を計算します。

# x軸の値を指定:(0と負の整数は計算できない)
x_vals = np.arange(start=-5.0, stop=5.0, step=0.002)
print(x_vals[:5])

# ディガンマ関数の計算
y_vals = digamma(x_vals)
print(y_vals[:5])

[-5.    -4.998 -4.996 -4.994 -4.992]
[          nan -498.28766529 -248.28144827 -164.94189785 -123.26901391]

　0のとき計算結果が-inf、負の整数のときNaNになります。

　グラフに表示しない値をマスクします。(マスクしておかないと不連続な点を繋いだグラフになります。)

# 閾値を指定
threshold = 10.0

# 閾値外の値をマスク
y_mask_vals = np.ma.masked_where((y_vals < -threshold) | (y_vals > threshold), y_vals)
print(y_mask_vals[:5])

[nan -- -- -- --]

　閾値thresholdを指定して、その範囲外の値の要素をnp.ma.masked_where()でマスクします。第1引数にマスクする条件、第2引数にマスクしない場合に返す値を指定します。この例では、$y < - \mathrm{threshold}$、$\mathrm{threshold} < y$の値をマスクします。複数の条件A, Bの一方でも該当する要素をマスクする場合は、(A) | (B)と書きます。

　ディガンマ関数のグラフを作成します。

# ディガンマ関数を作図
plt.figure(figsize=(12, 9)) # 図の設定
plt.plot(x_vals, y_mask_vals, color='orange') # 折れ線グラフ
plt.xlabel('x') # x軸ラベル
plt.ylabel('$\psi(x)$') # y軸ラベル
plt.suptitle('Digamma Function', fontsize=20) # 全体のタイトル
plt.xticks(ticks=np.arange(start=np.floor(x_vals.min()), stop=np.ceil(x_vals.max())+1.0)) # x軸目盛
plt.ylim(ymin=-threshold, ymax=threshold) # y軸の表示範囲
plt.grid() # グリッド線
plt.show() # 描画

ディガンマ関数のグラフ

　x軸の値が負の整数と0の点において、不連続なグラフになっているのが分かります。(マスクせずに線が繋がらないようにする引数の設定とかないですか？)

参考文献

• 岩田具治『トピックモデル』(機械学習プロフェッショナルシリーズ)講談社,2015年.

おわりに

　トリガンマ関数もあるみたいだけど、まだお目にかかったことがないのでやりません。

【次の内容】