【数学｜確率統計】任意の分布で乱数生成：逆関数法

概要

本記事は 任意の確率分布に基づく乱数生成方法 に関する記事です．

任意の確率分布に基づく乱数生成法として様々な方法が提案されている．

• 逆関数法(inversion method)
• 棄却法(reject method)
• 合成法(composition method)
• 合成棄却法(composition-rejection method)

本記事では 逆関数法 について説明し，具体例をいくつか示す．

---

逆関数法

確率変数$Y$が確率分布$f(y)$を持つとき,

X = \int_{-\infty}^{Y}dy f(y) \qquad (1)

の関係で定義される確率変数$X$は区間$[0, 1]$で一様分布する．
したがって，一様乱数$\{X_i\}$を発生させて，

\int_{-\infty}^{Y_i}dy f(y) = X_i

を満たす$\{ Y_i \}$を求めると，数列$\{ Y_i \}$は与えられた確率分布$f(y)$に従う乱数系列となる．

---

証明

式(1)より，$X$は$Y$の単調増加関数である．したがって，$Y$は$X$の1価関数で$Y = \phi(X)$と書くことができる．
したがって，

X = \phi^{-1}(Y)

が式(1)となる．この時，$Y$の分布は

\begin{align}
F_Y(y) &= P[Y < y] \\
    &= P[\phi(X) < y] \\
    &= P[X < \phi^{-1}(y)] \\
    &= \int_{0}^{\phi^{-1}(y)} dx f_{X}(x) \\
    &= \int_{0}^{\phi^{-1}(y)} dx \\
    &= \phi^{-1}(y)
\end{align}

となる．ここで，$f_{X}(x)$は$x$の分布であり，$X$は一様分布であるため$f_{X}(x) = 1$となる．
したがって、

F_{Y} = \int_{-\infty}^{y} d\eta f(\eta)

となる．
したがって，確率変数$Y$の確率分布は与えられた$f(y)$である．　[Q.E.D.]

---

例｜コーシー分布

コーシー分布(Cauchy distribution)は強制共鳴を記述する微分方程式の解や分光学のスペクトルの形状を記述するために用いられる分布である．

f(x; x_0, \gamma) = \frac{1}{\pi} \frac{\gamma}{(x - x_0)^2 + \gamma^2}

ここで，$x_0$は最頻値を与える位置母数，$\gamma$は半値半幅を与える尺度母数である．
式(1)を用いると，

\begin{align}
    X &= \int_{-\infty}^{Y} dy f(y; y_0, \gamma) \\
        &= \frac{1}{\pi} \arctan \left( \frac{Y-y_0}{\gamma} \right) + \frac{1}{2} \qquad (2)
\end{align}

となる．
したがって，確率変数Yについて解くと

\begin{align}
    Y = y_0 + \gamma \tan \left\{ \pi \left( X - \frac{1}{2} \right) \right\}
\end{align}

となる．
ここで，確率変数$X$は一様分布$U(0, 1)$に従う．

---

プログラム

以下のプログラムは$f(x; y_0 = 2, \gamma = 2)$に基づいた乱数を100000個生成する．

InversionMethod_cauchy.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# Import library
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats

# Setting parameter
np.random.seed(0)
N = 100000

nbins = 300000
y0 = 2
gamma = 2

# Create uniform random number
U = np.random.uniform(0.0, 1.0, N)

# Transform inversion
Y = y0 + gamma * np.tan( (U - 1/2) * np.pi )

# Plot result
plt.figure()
plt.hist(Y, nbins, density = True)

rv = scipy.stats.cauchy(loc = y0, scale = gamma)
x = np.linspace(rv.ppf(0.05), rv.ppf(0.95), 1000)
plt.plot(x, rv.pdf(x), 'r-', lw = 2)
plt.xlim(rv.ppf(0.05), rv.ppf(0.95))

plt.show()

---

結果

逆関数法を用いて生成した乱数のヒストグラムと確率分布$f(x; x_0, \gamma)$をプロットした結果を示す．
結果より，生成した乱数に基づくヒストグラムが確率分布に沿っていることが分かる．

---

あとがき

乱数生成の1つのである逆関数法を実装し，プラグラムを作成，乱数生成と確率分布をプロットし比較を行った．
次回は他の手法による乱数生成方法を実装してみることとする．

---

参考文献

1. 津田孝夫, 「モンテカルロ法とシミュレーション：電子計算機の確率論的応用」　, 培風館, 1995
2. yumi ito, "1. Pythonで学ぶ統計学　2. 確率分布[scipy.stats徹底理解]", https://qiita.com/y_itoh/items/c388ff82360906240daf, 2020, (サイトアクセス日：2024.5.31)
3. よし まる, "【Python】サンプリングアルゴリズム - 逆関数法", https://qiita.com/yoneXyone/items/9d7eb224ff19d6ea3d2e#%E4%BE%8B1%E6%8C%87%E6%95%B0%E5%88%86%E5%B8%83, 2024, (サイトアクセス日:2024.5.31)