この記事では、サンプリング手法のひとつである逆関数法について、その基本原理と実装例を紹介します。最近注目されている画像生成AIなどの分野でサンプリングの話題を目にする機会が増えていますが、実際の実装となると意外と複雑です。そこで、基本となる逆関数法に焦点をあて、具体例を通して理解を深めていきます。

※この記事作成にあたっては、こちらの記事を大変参考にさせて頂きました。

はじめに

「サンプリング」という言葉は、確率分布からサンプルを抽出する処理として広く使われています。私もかつては「母集団から確率に従ってサンプルを取り出す」というイメージだけで理解していたのですが、実際には多様な手法が存在し、画像生成AIの実装では拡散モデルといった複雑なサンプリング手法が採用されていることに気づかされました。
そこで今回は、基本中の基本となるサンプリング手法「逆関数法」について、基礎から実装例までを整理してみました。

逆関数法の基本原理

サンプリングの流れ

逆関数法は、次の3つのステップでサンプリングを実現します。

1.一様分布からの乱数生成
入力領域 [0, 1] の一様分布から乱数を生成します。

2.累積分布関数（CDF）の活用
生成した一様乱数は、CDF上の1点に対応付けられます。

3.逆関数による変換
CDFの逆関数に乱数を代入することで、目的の確率分布に従うサンプルを得ます。

この手法により、一様乱数から任意の確率分布に変換することが可能となります。

自作関数からのサンプリング

概要

サンプリングのイメージを掴むために、自作した関数を作ってサンプリングします。
具体的には、サイコロの目をサンプリングするシミュレーションを想定し、以下の手順でサンプリングします。

領域[0,1]を複数の領域に分割し、各領域の出力値を定義する
[0,1]の一様分布から乱数値を生成する
乱数値が入った領域の出力値を計算する

元々、入力は乱数つまり確率変数だったので、出力も確率変数となります。

今、出目の出現確率が等しくなるサイコロの目のシミュレーションを考えてみます。入出力のペアとして$\{(0〜1/6,1),(1/6〜2/6,2),(2/6〜3/6,3),(3/6〜4/6,4),(4/6〜5/6,5),(5/6〜6/6,6)\}$を考え、[0,1]の一様分布から乱数値を入力すると、一様分布の出力値となり、シミュレーションできそうですね。

では、某アニメのように4と5と6の目しかない出ないサイコロのシミュレーションはどうでしょうか？入出力のペアとして入出力のペアとして$\{(0〜2/6,4),(2/6〜4/6,5),(4/6〜6/6,6)\}$と考え、[0,1]の一様分布から乱数値を入力すると、シミュレーションできそうですね。当然先ほどの確率分布とは異なります。

このように定義域によって一様乱数から、任意の確率分布に変換することが可能となります。この時の入力の定義域は累積分布関数の出力値の定義域と等しくなります。そのため、逆関数法によってサンプリングしていると言えます。

実装：サイコロの目をサンプリングする

ライブラリのインポート

inverse_function.py

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats

サンプリング

inverse_function.py

# 乱数値が領域内か否かを判定する関数。
def section(u,a,b):
    if (a<=u)&(u<b):
        return 1
    else:
        return 0

# パラメータの設定
np.random.seed(0)
N = 100000
N_dice = 6

# 確率分布関数を rv に格納
# 範囲[0.0, 1.0]の一様分布からの確率変数
rv = uniform(loc = 0.0, scale = 1.0)
binarys = [[section(u,num/N_dice,(num+1)/N_dice) for num in range(0,N_dice)] for u in rv.rvs(size = N)]

結果の確認

inverse_function.py

# 出目の出現頻度を確認
freq = np.sum(binarys,axis=0)/N
print(freq)

確率分布関数からのサンプリング

既存の確率分布関数を利用する場合は、定義域は利用する関数の形に依存します。
指数分布関数をもとに確認していきます。

f(x|\lambda)=\lambda e^{-\lambda x}

この関数の累積分布関数は以下のようになります。

F(x|\lambda)=1-\lambda e^{-\lambda x}

もちろん$F(x|\lambda)=[0,1]$はです。
ここで[0,1]の一様乱数により値$u$が得られたとします。
この時以下の式が成り立ちます。

u=1-\lambda e^{-\lambda x}

この式を$x$の右辺が等式になるよるように変換すると、左辺の指揮は一様乱数から、指数分布関数に従う変数$x$が生成されるという意味になります。

x=F^{-1}(u)=-\frac{1}{\lambda} \log(1-u)

この一連の手続きが逆関数法であり、一様分布から既知の確率分布関数に従う変数を得ることができます。

実装例：指数分布関数からのサンプリング

ライブラリのインポート

inverse_function.py

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
# シードの固定
np.random.seed(0)

サンプリング

inverse_function.py

nbins = 50
# 指数分布のパラメータ(scale = 1/lambda)
scale = 1.0
N = 100000

#一様乱数を生成
U = scipy.stats.uniform(loc = 0.0, scale = 1.0).rvs(size = N)

# 指数分布の累積分布関数の逆関数を用いて変換
X1 = - scale * np.log(1-U)

結果の確認

inverse_function_draw.py

# 変換した指数分布の乱数と理想的な確率密度関数(PDF)を描画
plt.figure()

# 変換した分布関数のプロット
plt.hist(X1, int(nbins * 10 /5), density = True)

# 指数分布の確率密度関数を定義
rv = scipy.stats.expon(scale = scale)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)

# 真の確率密度関数(pdf)をプロット
plt.plot(x, rv.pdf(x), 'r-', lw = 2)

# 上端と下端を設定
plt.xlim(rv.ppf(0.01),rv.ppf(0.99))
plt.show()

おわりに

最後まで読んでいただきありがとうございました。

本記事を執筆する動機は、拡散モデルに向かって、一通りのサンプリング法の知識をまとめていくことにあります。逆関数法は、基礎的なサンプリング手法ですが、サンプリング＝利用しやい乱数から欲しい確率分布に従う乱数を生成すること、の理解を深めるにあたり良い題材になると、執筆して感じました。

図などを導入し、もう少し視覚的に読みやすくする工夫は今後の課題とします。

知識至らずで曖昧な点があると思います。ご指摘等ありましたら、編集いたしますので、コメントお待ちしております。