はじめに

　機械学習や統計学で登場する各種の確率分布について、「計算式の導出・計算のスクラッチ実装・計算過程や結果の可視化」などの「数式・プログラム・図」を用いた解説により、様々な角度から理解を目指すシリーズです。

　この記事では、ガウス分布の確率密度などの計算についてPythonを使って確認します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
1次元ガウス分布の計算
- 確率密度の計算
- 統計量の計算
  - スクラッチ実装による計算
  - クラスによる計算
参考文献
おわりに

1次元ガウス分布の計算

　1次元ガウス分布(Gaussian distribution・正規分布・Normal distribution)に関する計算を実装します。この記事では、Pythonを利用して計算します。
　ガウス分布については「1次元ガウス分布の定義式 - からっぽのしょこ」、Rを利用する場合は「【R】1次元ガウス分布の計算 - からっぽのしょこ」を参照してください。

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

# ライブラリを読込
import numpy as np
import scipy as sp

　また、スクラッチ実装で利用するクラスを読み込みます。

# クラスを読込
from scipy.stats import norm # ガウス分布

　この記事では(コードがごちゃごちゃしないように)、SciPyライブラリからクラスなどを個別に読み込んで利用します。

確率密度の計算

　まずは、1次元ガウス分布に従う確率密度(density)を計算します。

パラメータの設定

　ガウス分布のパラメータ $\mu, \sigma$ を設定します。

# 平均パラメータを指定
mu = 2.0
print(mu)

# 標準偏差パラメータを指定
sigma = 0.5
print(sigma)

2.0
0.5

　平均(実数) $\mu$ 、標準偏差(正の値) $\sigma \gt 0$ を指定します。

　また、計算式によって必要になるパラメータ $\lambda$ を設定します。

　標準偏差(分散)から精度を作成します。

# 精度パラメータに変換
lmd = 1.0 / sigma**2
print(lmd)

4.0

　標準偏差 $\sigma$ (分散 $\sigma^2$ )を用いて、精度 $\lambda = \frac{1}{\sigma^2}$ を計算します。

　または、精度から標準偏差を作成します。

# 精度パラメータを指定
lmd = 4.0
print(lmd)

# 標準偏差パラメータに変換
sigma = 1.0 / np.sqrt(lmd)
print(sigma)

4.0
0.5

　精度(正の値) $\lambda \gt 0$ を指定して、標準偏差 $\sigma = \frac{1}{\sqrt{\lambda}}$ を計算します。

　確率変数の実現値 $x$ を設定します。

# 確率変数の値を指定
x = 1.5
print(x)

1.5

　 $x$ を指定します。

　続いて、設定した値に従う確率密度を計算していきます。

スクラッチ実装による計算

　標準偏差・精度のパラメータに関してそれぞれ実装します。

標準偏差パラメータの場合

　定義式を実装して、確率密度を計算します。

# 定義式により確率密度を計算
dens  = 1.0 / np.sqrt(2.0 * np.pi * sigma**2)
dens *= np.exp(-0.5 * (x - mu)**2 / sigma**2)
print(dens)

0.48394144903828673

　1次元ガウス分布は、次の式で定義されます。

$\displaystyle \begin{aligned} \mathrm{C}_{\mathcal{N}} &= \frac{1}{\sqrt{2 \pi \sigma^2}} \\ \mathcal{N}(x \mid \mu, \sigma^2) &= \mathrm{C}_{\mathcal{N}} \exp \Bigl( - \frac{1}{2 \sigma^2} (x - \mu)^2 \Bigr) \end{aligned}$

　ここで、 $\mathrm{C}_{\mathcal{N}}$ はガウス分布の正規化項です。
　円周率 $\pi$ は np.pi で扱えます。
　ネイピア数 $e$ の指数関数 $e^x = \exp(x)$ は、np.exp() で計算できます。

　対数をとった定義式を実装して、確率密度を計算します。

# 定義式により確率密度を計算
log_dens  = -0.5 * np.log(2.0 * np.pi) - np.log(sigma)
log_dens -= 0.5 * (x - mu)**2 / sigma**2
dens      = np.exp(log_dens)
print(log_dens)
print(dens)

-0.7257913526447274
0.48394144903828673

　対数をとった1次元ガウス分布は、次の式になります。

$\displaystyle \begin{aligned} \log \mathrm{C}_{\mathcal{N}} &= - \frac{1}{2} \log (2 \pi) - \log \sigma \\ \log \mathcal{N}(x \mid \mu, \sigma^2) &= \log \mathrm{C}_{\mathcal{N}} - \frac{1}{2 \sigma^2} (x - \mu)^2 \end{aligned}$

　指数と対数の関係 $\exp(\log x) = x$ より、計算結果の指数をとると確率密度が得られます。

$\displaystyle \mathcal{N}(x \mid \mu, \sigma^2) = \exp \Bigl( \log \mathcal{N}(x \mid \mu, \sigma^2) \Bigr)$

　自然対数 $\log x$ は np.log() で計算できます。

精度パラメータの場合

　続いて、精度 $\lambda$ を用いて、確率密度を計算します。

# 定義式により確率密度を計算
dens  = np.sqrt(0.5 * lmd / np.pi)
dens *= np.exp(-0.5 * lmd * (x - mu)**2)
print(dens)

0.48394144903828673

　パラメータに精度 $\lambda$ を用いると、次の式になります。

$\displaystyle \begin{aligned} \mathrm{C}_{\mathcal{N}} &= \sqrt{\frac{\lambda}{2 \pi}} \\ \mathcal{N}(x \mid \mu, \lambda^{-1}) &= \mathrm{C}_{\mathcal{N}} \exp \Bigl( - \frac{\lambda}{2} (x - \mu)^2 \Bigr) \end{aligned}$

　対数をとった定義式を実装して、確率密度を計算します。

# 定義式により確率密度を計算
log_dens  = 0.5 * (np.log(lmd) - np.log(2.0 * np.pi))
log_dens -= 0.5 * lmd * (x - mu)**2
dens      = np.exp(log_dens)
print(log_dens)
print(dens)

-0.7257913526447274
0.48394144903828673

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

$\displaystyle \begin{aligned} \log \mathrm{C}_{\mathcal{N}} &= \frac{1}{2} \Bigl( \log \lambda - \log (2 \pi) \Bigr) \\ \log \mathcal{N}(x \mid \mu, \lambda^{-1}) &= \log \mathrm{C}_{\mathcal{N}} - \frac{\lambda}{2} (x - \mu)^2 \end{aligned}$

クラスによる計算

　標準偏差・精度のパラメータ、また標準化した変数に関してそれぞれ実装します。

標準偏差パラメータの場合

　scipy のクラスメソッドを使って、確率密度を計算します。

# クラスにより確率密度を計算
dens = norm.pdf(x=x, loc=mu, scale=sigma)
print(dens)

0.48394144903828673

　1次元ガウス分布の確率密度は、sp.stats.norm.pdf() で計算できます。確率変数の引数 x に $x$ 、平均の引数 loc に $\mu$ 、標準偏差の引数 scale に $\sigma$ を指定します。

　logpdf() を使うと、対数をとった確率密度を計算します。

# クラスにより確率密度を計算
log_dens = norm.logpdf(x=x, loc=mu, scale=sigma)
dens     = np.exp(log_dens)
print(log_dens)
print(dens)

-0.7257913526447274
0.48394144903828673

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

精度パラメータの場合

　続いて、精度 $\lambda$ を用いて、確率密度を計算します。

# クラスにより確率密度を計算
dens = norm.pdf(x=x, loc=mu, scale=1.0/np.sqrt(lmd))
print(dens)

0.48394144903828673

　標準偏差 $\sigma = \frac{1}{\sqrt{\lambda}}$ に変換して scale 引数に指定します。

標準化変数の場合

　標準化した変数 $z$ を用いて、確率密度を計算します。

# 変数を標準化
z = (x - mu) / sigma
print(z)

# 関数により確率密度を計算
dens = norm.pdf(x=z, loc=0.0, scale=1.0) / sigma
print(dens)

-1.0
0.48394144903828673

　パラメータ $\mu, \sigma$ を用いて、変数 $x$ を標準化 $z = \frac{x - \mu}{\sigma}$ します。
　標準化した変数 $z$ を用いると、次の式になります。

$\displaystyle \mathcal{N}(x \mid \mu, \sigma^2) = \frac{1}{\sigma} \mathcal{N}(z \mid 0, 1^2)$

　平均 $0$ 、標準偏差 $1$ (分散 $1$ )のガウス分布 $\mathcal{N}(x \mid \mu = 0, \sigma^2 = 1)$ を標準ガウス分布(標準正規分布) $\mathcal{N}(x \mid 0, 1)$ と呼びます。
　式変形については「分布の定義式」を参照してください。

　ここまでで、1次元ガウス分布の確率の計算を確認しました。次からは、分布の特徴などを表す数値の計算を確認していきます。

スポンサードリンク

統計量の計算

　次は、1次元ガウス分布の統計量(statistics)を計算します。
　計算式については「準備中」を参照してください。

スクラッチ実装による計算

　計算式を実装して、期待値を計算します。

# 計算式により期待値を計算
E_x = mu
print(E_x)

2.0

　1次元ガウス分布の期待値は、定義より、平均パラメータ $\mu$ です。

$\displaystyle \mathbb{E}[x] = \mu$

　分散を計算します。

# 計算式により分散を計算
V_x = sigma**2
print(V_x)

0.25

　分散は、定義より、標準偏差パラメータ $\sigma$ の2乗です。

$\displaystyle \mathbb{V}[x] = \sigma^2$

　最頻値を計算します。

# 計算式により最頻値を計算
mode_x = mu
print(mode_x)

2.0

　最頻値(モード)は、平均パラメータ $\mu$ です。

$\displaystyle \mathrm{mode}[x] = \mu$

クラスによる計算

　scipy のクラスメソッドを使って、期待値を計算します。

# クラスにより期待値を計算
E_x = norm.mean(loc=mu)
print(E_x)

2.0

　1次元ガウス分布の期待値は、sp.stats.norm.mean() で計算できます。確率密度メソッド pdf() と同様に引数を指定します。

　分散を計算します。

# 関数により分散を計算
V_x = norm.var(scale=sigma)
print(V_x)

0.25

　分散は、分散メソッド var() で計算できます。

　この記事では、1次元のガウス分布に関する計算を確認しました。次の記事では、確率分布のグラフを作成します。

参考文献

機械学習スタートアップシリーズ　ベイズ推論による機械学習入門

作者:須山敦志
講談社

Amazon

おわりに

2025.01.06：「1次元ガウス分布の作図」から記事を独立して、加筆修正しました。

　定義上仕方がない訳ですが、他の分布の記事と比べるとモーメントの類に解説記事としての面白みが足りないですね。せめてもの抵抗の意味も込めて(いや書いてから更新までに間が空いてしまったので当初の意図は忘れてしまいましたが)、標準化した計算などを追加しましたが、標準正規分布の計算しかできない状況が思い当たらないので必要性は感じませんね。

　最後に、先日公開されたukkaの新曲のライブ映像をどうぞ♪

　先日グループの解散が発表されまして、、、新メンバーもすっかり馴染んで各メンバーの対外的な活動も実ってきて(るようににわかオタクな私にも見えて)、今後のグループや個人の活躍に更なる期待がされている頃合いだろうに、どうして、、、むしろこれからもっと活動を拡げていくように感じていたのに、、、
　個人的にも「Aonity」で更にukkaにハマったと感じる新曲だったのに、本当にどうしてなんだ、、、

　せめてここまで辿り着いた皆さん何曲か聴いてみてください、もう半年足らずの活動期間なんです。お願いします。

【次の内容】

からっぽのしょこ

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

【Python】1次元ガウス分布の計算

はじめに

1次元ガウス分布の計算

確率密度の計算

パラメータの設定

スクラッチ実装による計算

標準偏差パラメータの場合

精度パラメータの場合

クラスによる計算

標準偏差パラメータの場合

精度パラメータの場合

標準化変数の場合

統計量の計算

スクラッチ実装による計算

クラスによる計算

参考文献

おわりに