はじめに

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

　この記事では、ガウス分布の乱数についてPythonを使って確認します。

【前の内容】

zenn.dev

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
1次元ガウス分布の乱数生成
- サンプリング
- サンプルサイズの影響
  - サンプルサイズと形状の関係
参考文献
おわりに

1次元ガウス分布の乱数生成

　1次元ガウス分布(Gaussian distribution・正規分布・Normal distribution)に従う乱数(random number)を生成して、グラフやアニメーションで確認します。この記事では、PythonのMatplotlibライブラリを利用して作図します。
　ガウス分布については「1次元ガウス分布の定義式 - からっぽのしょこ」、グラフ作成については「【Python】1次元ガウス分布の作図 - からっぽのしょこ」、Rを利用する場合は「【R】1次元ガウス分布の乱数生成 - からっぽのしょこ」を参照してください。

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

# ライブラリを読込
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

サンプリング

　まずは、1次元ガウス分布に従う乱数を作成します。

パラメータの設定

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

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

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

0.0
2.5

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

　続いて、設定した値に従う乱数を生成して、グラフを作成していきます。

乱数の生成

　サンプルサイズ $N$ を指定して、ガウス分布の乱数 $x_n$ を生成します。

# サンプルサイズを指定
N = 1000

# ガウス分布に従う乱数を生成
x_n = np.random.normal(loc=mu, scale=sigma, size=N)
print(x_n[:5])

[ 3.59332552 -2.00039499  5.17217649  0.21193451  2.15583667]

　サンプルサイズ(データ数) $N$ を指定します。
　 $N$ 個のサンプル(乱数)の集合 $\mathbf{x} = \{x_1, \cdots, x_N\}$ を作成します。
　1次元ガウス分布の乱数は、np.random.normal() で生成できます。平均の引数 loc に $\mu$ 、標準偏差の引数 scale に $\sigma$ 、サンプルサイズの引数 size に $N$ を指定します。

　以上で、乱数を作成できました。続いて、乱数を可視化していきます。

乱数の作図

　確率変数 $x$ の集計範囲を設定します。

# 階級値の範囲を設定
#x_min = -10.0
#x_max = 10.0
u = 5.0
x_size = np.max(np.abs(x_n - mu)) # 基準値を指定
x_min  = mu - x_size
x_max  = mu + x_size
x_min  = np.floor(x_min /u)*u # u単位で切り下げ
x_max  = np.ceil(x_max /u)*u  # u単位で切り上げ
print(x_min, x_max)

-10.0 10.0

　集計や作図に用いるため、階級値の最小値・最大値を作成します。
　この例では、生成したサンプル x_n の最小値・最大値を使って、範囲を設定しています。

　ガウス分布の乱数のヒストグラムを作成します。

# ラベル用の文字列を作成
param_lbl = f'$N = {N}, \\mu = {mu}, \\sigma = {sigma}$'

# サンプルの度数を作図
plt.figure(figsize=(9, 6), dpi=100, facecolor='white')
plt.hist(
    x=x_n, 
    bins=30, 
    color='#00A968'
) # 度数
plt.xlabel('$x$') # x軸ラベル
plt.ylabel('frequency') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('Gaussian distribution', fontsize=20) # 図タイトル
plt.grid() # グリッド線
plt.xlim(xmin=x_min, xmax=x_max) # x軸の表示範囲
plt.show()

　ヒストグラムは、plt.hist() で描画できます。データの引数 x にサンプルデータ、バーの数の引数 bins に(自動で設定される集計範囲における)バーの数を指定します。

　上の例では、階級数(バーの数)を指定して、階級幅や階級値は自動で設定されます。
　次の例では、階級幅や階級値を指定して、任意のバーを作図します。

　階級数から階級幅を作成します。

# 階級数を指定
bin_num = 20

# 階級幅を計算
bin_size = (x_max - x_min) / bin_num
print(bin_size)

1.0

　階級数 bin_num を用いて、階級幅 bin_size を計算します。
　階級値の最小値を $x_{\min}$ 、最大値を $x_{\max}$ として、階級幅(1区間のサイズ) $w$ は、階級値の範囲 $x_{\max} - x_{\min}$ を階級数(区間の数) $m$ で割った値 $w = \frac{x_{\max} - x_{\min}}{m}$ により求まります。

　または、階級幅から階級数を作成します。

# 階級幅を指定
bin_size = 1.0

# 階級数を計算
bin_num = ((x_max - x_min) // bin_size).astype(dtype='int')
print(bin_num)

　階級幅 bin_size を用いて、階級数 bin_num を計算します。
　階級数(区間の数) $m$ は、階級値の範囲 $x_{\max} - x_{\min}$ を階級幅(1区間のサイズ) $w$ で割った整数 $m = \lfloor \frac{x_{\max} - x_{\min}}{w} \rfloor$ により求まります。
　 $\frac{x_{\max} - x_{\min}}{w}$ が割り切れない(整数にならない)場合は、 $w$ を再設定する必要があります。
　作図の都合上、階級数(バーの数)を整数型にしておきます。

　階級値と境界値を作成します。

# 境界値の範囲を設定
bin_min = x_min - 0.5*bin_size
bin_max = x_max + 0.5*bin_size

# 階級値を作成
center_vec = np.arange(start=x_min, stop=x_max+bin_size, step=bin_size)
print(center_vec[:5])
print(center_vec.shape)

# 境界値を作成
bin_vec = np.arange(start=bin_min, stop=bin_max+0.5*bin_size, step=bin_size) # (x+0.5*wはstop引数が未満のための対策用)
#bin_vec = np.hstack([center_vec, center_vec[-1]+bin_size]) - 0.5*bin_size
print(bin_vec[:5])
print(bin_vec.shape)

[-10.  -9.  -8.  -7.  -6.]
(21,)
[-10.5  -9.5  -8.5  -7.5  -6.5]
(22,)

　階級値の最小値 $x_{\min}$ から最大値 $x_{\max}$ までを階級幅 $w$ の間隔で分割して、各区間の階級値(境界値の中央値)を作成します。
　境界値(区間全体)の最小値 $x_{\min} - \frac{w}{2}$ から最大値 $x_{\max} + \frac{w}{2}$ までを階級幅 $w$ の間隔で分割して、各区間の境界値(階級値の中央値)を作成します。または、各境界値から階級幅の半分 $\frac{w}{2}$ をスライドした値を作成します。

　階級を指定して、ヒストグラムを作成します。

# サンプルの度数を作図
plt.figure(figsize=(9, 6), dpi=100, facecolor='white')
plt.hist(
    x=x_n, 
    bins=bin_num+1, range=(bin_min, bin_max), 
    color='#00A968'
) # 度数
plt.xticks(ticks=center_vec) # 階級値
plt.xlabel('$x$') # x軸ラベル
plt.ylabel('frequency') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('Gaussian distribution', fontsize=20) # 図タイトル
plt.grid() # グリッド線
plt.show()

　バーの数の引数 bins に階級数、バーの範囲の引数 range に境界値(区間全体)の最小値・最大値を指定します。バーの数や幅は range 引数の設定により決まります。

　ガウス分布に従う確率密度を計算します。装飾用の処理です。

# x軸の値を作成
x_vec = np.linspace(start=bin_min, stop=bin_max, num=1001)

# ガウス分布の確率密度を計算
gen_dens_vec = norm.pdf(x=x_vec, loc=mu, scale=sigma)

　詳しくは「分布の作図」を参照してください。

　密度のグラフを作成します。

# サンプルの度数を作図
plt.figure(figsize=(9, 6), dpi=100, facecolor='white')
plt.hist(
    x=x_n, 
    bins=bin_num+1, range=(bin_min, bin_max), density=True, 
    color='#00A968', alpha=0.5, 
    label='random number'
) # 密度
plt.plot(
    x_vec, gen_dens_vec, 
    color='green', linewidth=1.0, linestyle='--', 
    label='generator'
) # 確率密度
plt.xticks(ticks=center_vec) # 階級値
plt.xlabel('$x$') # x軸ラベル
plt.ylabel('density') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('Gaussian distribution', fontsize=20) # 図タイトル
plt.legend(title='distribution') # 凡例
plt.grid() # グリッド線
plt.show()

　形状の比較用に、ガウス分布に従う確率密度(サンプルの生成分布)を破線で示します。

　density=True を指定すると、縦軸を密度に変換したヒストグラムを描画します。

　ここまでで、作図関数内の集計処理によって、グラフを作成できました。続いては、データフレーム上での集計処理によって、集計結果と共に同様のグラフを作成します。

乱数の集計

　サンプル集合 $\mathbf{x}$ を集計します。

# 度数を集計
freq_vec, bin_vec = np.histogram(a=x_n, bins=bin_num+1, range=(bin_min, bin_max))
print(bin_vec[:10])
print(freq_vec[:10])

# 相対度数を計算
relfreq_vec = freq_vec / N
print(relfreq_vec[:10])

# 密度を計算
sample_dens_vec = freq_vec / (bin_size * N)
print(sample_dens_vec[:10])

# 階級値を作成
center_vec = bin_vec[:-1] + 0.5*bin_size
print(center_vec[:10])

[-10.5  -9.5  -8.5  -7.5  -6.5  -5.5  -4.5  -3.5  -2.5  -1.5]
[  0   0   0   1  11  20  44  84 109 145]
[0.    0.    0.    0.001 0.011 0.02  0.044 0.084 0.109 0.145]
[0.    0.    0.    0.001 0.011 0.02  0.044 0.084 0.109 0.145]
[-10.  -9.  -8.  -7.  -6.  -5.  -4.  -3.  -2.  -1.]

　サンプルの度数は、np.histogram() でカウントできます。データの引数 a にサンプルデータ、バーの数の引数 bins に階級数、バーの範囲の引数 range に境界値の最小値・最大値を指定します。
　度数を $N_x$ ( 階級値が $1$ の範囲の度数を $N_1$ で表す)として、サンプルサイズ $N$ で割ると相対度数 $\frac{N_x}{N}$ 、更に階級幅 $w$ で割ると密度 $\frac{N_x}{w N}$ が求まります。

　または、density=True を指定すると、密度を出力します。

# 密度を集計
smp_dens_vec, bin_vec = np.histogram(a=x_n, bins=bin_num+1, range=(bin_min, bin_max), density=True)
print(bin_vec[:10])
print(smp_dens_vec[:10])

[-10.5  -9.5  -8.5  -7.5  -6.5  -5.5  -4.5  -3.5  -2.5  -1.5]
[0.    0.    0.    0.001 0.011 0.02  0.044 0.084 0.109 0.145]

　ガウス分布の乱数のヒストグラムを作成します。

# サンプルの度数を作図
plt.figure(figsize=(9, 6), dpi=100, facecolor='white')
plt.bar(
    x=center_vec, height=freq_vec, 
    width=bin_size, align='center', 
    color='#00A968'
) # 度数
#plt.xticks(ticks=center_vec) # 階級値
plt.xlabel('$x$') # x軸ラベル
plt.ylabel('frequency') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('Gaussian distribution', fontsize=20) # 図タイトル
plt.grid() # グリッド線
plt.show()

　棒グラフの描画関数 plt.bar() でヒストグラムを作成します。横軸の引数 x に階級値、縦軸の引数 height に度数を指定します。また、(棒グラフではなく)ヒストグラムの体裁にするには、バーの幅の引数 width に階級幅を指定します。

　密度のグラフを作成します。

# サンプルの密度を作図
plt.figure(figsize=(9, 6), dpi=100, facecolor='white')
plt.bar(
    x=center_vec, height=smp_dens_vec, width=bin_size, 
    color='#00A968', alpha=0.5, 
    label='random number'
) # 密度
plt.plot(
    x_vec, gen_dens_vec, 
    color='green', linewidth=1.0, linestyle='--', 
    label='generator'
) # 確率密度
#plt.xticks(ticks=center_vec) # 階級値
plt.xlabel('$x$') # x軸ラベル
plt.ylabel('density') # y軸ラベル
plt.title(param_lbl, loc='left') # タイトル
plt.suptitle('Gaussian distribution', fontsize=20) # 図タイトル
plt.legend(title='distribution') # 凡例
plt.grid() # グリッド線
plt.show()

　height 引数に密度(ヒストグラム全体に対するバーの面積の割合)を指定します。

　第2軸(右側)に度数目盛を表示するグラフを作成します。

# 密度軸の範囲を指定
u = 0.05
dens_max = np.max(freq_vec / (bin_size * N))
dens_max = np.ceil(dens_max /u)*u # u単位で切り上げ

# サンプルの密度を作図
fig, ax = plt.subplots(figsize=(9, 6), dpi=100, facecolor='white')
fig.suptitle('Gaussian distribution', fontsize=20)

# 第1軸の図
ax.bar(
    x=center_vec, height=smp_dens_vec, width=bin_size, 
    color='#00A968', alpha=0.5, 
    label='random number'
) # 密度
ax.plot(
    x_vec, gen_dens_vec, 
    color='green', linewidth=1.0, linestyle='--', 
    label='generator'
) # 確率密度
#ax.set_xticks(ticks=center_vec) # 階級値
ax.set_xlabel('$x$')
ax.set_ylabel('density') # 1軸ラベル
ax.set_title(param_lbl, loc='left')
ax.legend(title='distribution') # 凡例
ax.grid() # グリッド線
ax.set_ylim(ymin=0.0, ymax=dens_max) # (目盛の共通化用)

# 2軸用のインスタンスを作成
ax2 = ax.twinx()

# 密度軸(1軸)目盛を取得
dens_vals = ax.get_yticks()

# 度数軸(2軸)目盛に変換
freq_vals = dens_vals * bin_size * N

# 第2軸の図
ax2.set_yticks(ticks=freq_vals, labels=[f'{y:.1f}' for y in freq_vals]) # 2軸目盛
ax2.set_ylabel('frequency') # 2軸ラベル
ax2.set_ylim(ymin=0.0, ymax=dens_max*bin_size*N) # (目盛の共通化用)

plt.show()

　(2軸の図は、pyplot.function() の記法では作成できないようなので、axes.function() の記法で作成します。)

　第2軸は、ax.twinx() を使って描画できます。x軸を共有するの描画用のインスタンスが出力して、第1軸と同様に軸の設定を追加します。
　2つの軸の目盛の位置を共通化するために、第1軸の目盛の値を ax.get_yticks() で取り出して、第2軸の目盛の値として ax.set_yticks() で設定します。この例では、第1軸は密度なので、階級幅とデータ数(サンプルサイズ)を掛けて度数に変換して第2軸に使います。
　軸の範囲を ax.set_ylim() で対応させます。

　以上で、基本的な乱数生成の処理を確認しました。

スポンサードリンク

サンプルサイズの影響

　次は、サンプルサイズを増やしてグラフの変化をアニメーションで確認します。
　作図コードについては「Probability-Distribution/code/gaussian/random_number.py at main · anemptyarchive/Probability-Distribution · GitHub」を参照してください。

サンプルサイズと形状の関係

　サンプルサイズ $N$ を大きくしたときのガウス分布の乱数のヒストグラム(度数)の変化をアニメーションにします。

　上の図は1サンプルずつ、下の図は複数サンプルずつ増やしたときの様子です。

　密度の変化をアニメーションにします。

　サンプルが増えるほど、生成分布の形状に近付くのが分かります。

　以上で、サンプルサイズの影響を確認しました。

　この記事では、1次元のガウス分布の乱数を生成しました。次の記事では、生成した乱数をガウス分布のパラメータとして利用します。

参考文献

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

作者:須山敦志
講談社

Amazon

おわりに

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

　初稿時と比較すると、第2軸を装飾できるようになりました。

　2025年1月5日は、超ときめき♡宣伝部の小泉遥香さんの25歳のお誕生日です。

　恵比寿方面で楽しんでいたらおはるさん&とき宣に辿り着きました、宣伝部員を名乗るにはまだ早い、ど新規です。よろしくお願いします。
　おはるさんのお顔がとにかくどタイプでして、SNSなどでぼんやり眺めていたらガハハと笑う姿にも惹かれ、歌う姿にも魅せられました。とき宣の楽曲についても、今どきのショート動画映えするキャッチーなカワイイ曲だけでなく、ライブで盛り上がりそうな格好いい曲もあり、ハマってる途中です♬
　とき宣としてだけでも10年以上もの期間を活動してくれていたおかげで巡り会えました、ありがとうございます。
　そして、お誕生日おめでとうございます！これからの活躍をリアルタイムで観られるのを楽しみにしています。

【次の内容】

からっぽのしょこ

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

【Python】1次元ガウス分布の乱数生成