任意の確率密度関数を持つ乱数の生成

はじめに

この記事では、Pythonを用いて任意の確率密度関数に従う乱数を生成する方法について解説します。特に、三角形内部に均一に乱数を配置する問題を例に、具体的なコードと数式を用いて丁寧に説明します。累積分布関数とその逆関数を用いることで、一様乱数からどのようにして任意の確率密度関数を持つ乱数を生成できるのかを理解することができます。

対象読者:

• 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$ が出現するようにします。これは、$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}$ の導出方法を説明します。

1. 累積分布関数の計算: 確率密度関数 $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$であることに注意すると、累積分布関数のグラフは以下のようになります。

2. 逆関数の計算: 累積分布関数 $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$ に変換する関数は、累積分布関数の逆関数となるのです。

つまり、累積分布関数の逆関数は、一様分布を目的の分布に「歪める」役割を果たします。