一様乱数から任意の分布を作る：逆関数法の仕組みと実装

はじめに

三角形の内部に点を一様にばらまきたいとき、頂点から2つの一様乱数で線形補間するだけでは、点が三角形の一部に偏ってしまいます。均一な散らばりを得るには、使う乱数の分布そのものを作り替えなければなりません。この三角形の問題を手がかりに、一様乱数から任意の確率密度関数に従う乱数を作る手順を追います。鍵になるのは累積分布関数（cumulative distribution function）とその逆関数です。

対象読者:

• Pythonの基本的なプログラミングスキルを持つ方
• 確率・統計の基礎知識（確率密度関数、累積分布関数など）を持つ方
• 乱数生成のアルゴリズムに興味がある方

記事のポイント:

• 一様乱数から任意の確率密度関数に従う乱数を作る一般的な手順を解説
• 累積分布関数とその逆関数が果たす役割を整理
• 三角形の内部に点を均一に配置するPythonコードを提示
• 変換式の導出を数式で確認

一様乱数から始める

乱数生成の基本は、一様乱数です。一様乱数とは、ある範囲内のすべての値が等しい確率で出現する乱数のことです。

確率密度関数で表すと、[0, 1) の範囲で一様乱数の確率密度関数 $P(x)$ は1になります。

$$P(x) = 1, \quad 0 \le x < 1$$

Pythonでは、randomモジュールを使って一様乱数を生成できます。

# [0, 1)の乱数を生成
import random
print(random.random())

定数倍と平行移動で、範囲の異なる一様乱数も作れます。

# [0, 2)の乱数を生成
print(random.random() * 2)

# [-2, 2)の乱数を生成
print(random.random() * 4 - 2)

三角形の内部に点を均一にばらまく

冒頭の三角形の問題に戻り、乱数の分布を作り替える手順を具体的に見ていきます。

線形補間だけでは偏る

三角形の頂点をp0, p1, p2とし、それぞれの座標を次のように定めます。

• p0 = (0, 0, 0)
• p1 = (1, 0, 0)
• p2 = (0, 1, 0)

まず、素朴に線形補間だけで点を打つコードです。

import matplotlib.pyplot as plt
import random
import numpy as np

p0 = np.array([0,0,0])
p1 = np.array([1,0,0])
p2 = np.array([0,1,0])

v01 = p1 - p0
v02 = p2 - p0

points = []
for i in range(400):
    t = random.random()
    u = random.random()
    points.append(p0 + t * v01 + (1 - t) * u * v02)
points = np.array(points)

plt.figure(figsize=(5,5))
plt.scatter(points[:, 0], points[:,1])
plt.show()

このコードでは2つの乱数 $t$ と $u$ を用いて三角形内部の点を計算しています。しかし、この方法では均一な分布にはなりません。$t$ が一様分布であるため、p0に近い部分と遠い部分で $v_{02}$ 方向のスパンが異なるにもかかわらず、同じ数の点がそのスパン上に分布してしまうためです。

tの分布を作り替える

偏りをなくすには、$t$ の分布そのものを変えます。$t$ の位置での縦方向のスパンの長さに比例した頻度で $t$ が出るようにすればよく、そのスパンの長さは $1-t$ ですから、$1-t$ に比例した確率で $t$ が現れるようにします。

これを確率密度関数で書き、積分値が1になるように正規化すると $P(x) = 2(1-x)$ です。

この分布に従う乱数は、次のように生成できます。

import matplotlib.pyplot as plt
import random
import numpy as np

p0 = np.array([0,0,0])
p1 = np.array([1,0,0])
p2 = np.array([0,1,0])

v01 = p1 - p0
v02 = p2 - p0

points = []
for i in range(400):
    t = random.random()
    t = 1 - np.sqrt(1-t) # 変換
    u = random.random()
    points.append(p0 + t * v01 + (1 - t) * u * v02)
points = np.array(points)

plt.figure(figsize=(5,5))
plt.scatter(points[:, 0], points[:,1])
plt.show()

$t = 1 - \sqrt{1-t}$ と変換した $t$ を使うと、点は三角形の内部に均一に散らばります。

変換式の導出

変換式 $t := 1-\sqrt{1-t}$ は、累積分布関数を求めてから逆関数を解く、という2段階で導けます。

まず、確率密度関数 $P(x) = 2(1-x)$ を積分して累積分布関数 $f(x)$ を求めます。

$$f(x) = \int P(x) dx = \int 2(1-x) dx = 2x - x^2 + C$$

ここで、$f(0) = 0$ より、積分定数 $C = 0$ です。$x$が0未満、または1以上の時、$P(x)=0$であることに注意すると、累積分布関数のグラフは以下のようになります。

次に、この累積分布関数 $f(x) = 2x - x^2$ の逆関数を求めます。$y = 2x - x^2$ を $x$ について解きます。

$$x^2 - 2x + y = 0$$

二次方程式の解の公式より、

$$x = \frac{-(-2) \pm \sqrt{(-2)^2 - 4 \cdot 1 \cdot y}}{2 \cdot 1} = 1 \pm \sqrt{1 - y}$$

$0 ≤ x ≤ 1$ の範囲を考慮すると、

$$f^{-1}(y) = 1 - \sqrt{1 - y}$$

となります。これが求める変換式です。

なぜ逆関数で目的の分布になるのか

累積分布関数の逆関数を一様乱数に通すと目的の分布が得られる理由を、離散的な例で確かめます。100個の要素を持つ一様数列 $a = [0.00, 0.01, 0.02, ..., 0.99]$ を考え、この数列の各要素に、ある単調増加関数を適用した変換後の数列 $b$ が、特定の確率密度関数に従うようにしたいとします。

この「ある関数」が、数列 $b$ の累積分布関数 $f(x)$ の逆関数 $f^{-1}(x)$ です。累積分布関数 $f(x)$ は、値 $x$ が数列 $b$ の中で先頭から何%の位置にあるかを返します。一様数列 $a$ の各要素は、まさにこの「先頭から何%か」に当たる値です。位置から値を引き当てる関数、つまり累積分布関数の逆関数を適用すれば、数列 $b$ が得られます。

累積分布関数の逆関数は、一様分布を目的の分布に「歪める」役割を果たします。この手順は逆関数法（inverse transform sampling）と呼ばれ、一様乱数さえ生成できればどんな分布にも適用できる一般的な方法です。

まとめ

目的の確率密度関数を積分して累積分布関数を求め、その逆関数を一様乱数に適用すれば、任意の分布に従う乱数を作れます。三角形の例では、密度 $P(x) = 2(1-x)$ から変換式 $1-\sqrt{1-x}$ が閉じた形で求まりました。

ただし、逆関数がいつでも解析的に解けるとは限りません。解けない場合は累積分布関数を数値的に逆引きするか、正規分布のように専用のアルゴリズムが整っている分布なら numpy.random などライブラリの実装をそのまま使うのが実際的です。自分の問題で密度を指定したくなったら、まずその累積分布関数が積分できるか、逆関数が閉じた形で解けるかを確かめるところから始めてください。