Pythonでの生成方法を書いていきますが、根底にある考え方は他の言語でも共通です。
手始めに、一様乱数から
確立密度関数で表現すると、[0, 1)において、P(x) = 1。
# [0, 1)の乱数を生成
import random
print(random.random())
一様乱数を線形にスケールすることで、異なる範囲の一様乱数を生成できます。
# [0, 2)の乱数を生成
print(random.random() * 2)
# [-2, 2)の乱数を生成
print(random.random() * 4 - 2)
本題の、任意の確率密度関数を持つ乱数の生成方法
よくある（と私が思う）のは、三角形の内部に乱数で点を打ちたい、というシチュエーションです。例えば、
import matplotlib.pyplot as plt
import random
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()
として、生成するのが挙げられると思います。このとき、t + (1-t)uは0~1の範囲だから、pointsは、三角形の内部に分布します。ただし、この分布は均一ではありません。実際に描画してみると下図のようになります。なぜなら、まずv01(右方向)をt倍した位置から、v02(上方向)方向に(1-t)u倍した分だけ移動する、という意味になり、tが一様分布だと、p0に近い部分と遠い部分ではv02方向のスパンが異なるにもかかわらず同じ数の点がそのスパン上に分布することが乱数の分布として期待されるからです。
さて、これは、tに一様乱数を使わずに、別の分布を使えば解決します。いま問題になっているのは、v01のスケールtが一様乱数になっているせいで不当にp1周辺の点の密度が高くなっているのが原因と見ることができます。言い換えると、小さいtが現れやすく、1に近いtが現れにくい分布に置き換えたいわけです。
では、それがどのような分布かというと、tの位置の縦方向のスパンの長さに比例した頻度でtが出現してほしい、具体的な例を出すと、t=0.25の出現確率がt=0の出現確率の3/4になっていてほしい、(この例で)一般的に言うと、1-tに比例した確率でtが現れてほしい、ということになります。確立密度関数にすると、(積分が１となる制約が入りますので、) P(x) = 2(1-x) となります。そのような分布の乱数は次のように生成できます。
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)
きれいに分布していますね！
$$t := 1-\sqrt{1-t} $$
の変換を施しています。
この変換式をどうやって作るか、というのが問題ですが、step1として、確立密度関数を積分し、累積分布関数を計算します。
$$ f(x) = \int P(x) dx = \int 2(1-x) dx = 2x – x^2 \ $$
注意点は、０未満の範囲と、1以上の範囲のP(x)は0なので、グラフの形は次のようになります。
step2として、これの逆関数を作ります。累積分布関数を、xについて解きます。
$$ f^{-1}(x) =
\begin{cases}
1-\sqrt{1-x} &(0<x<1) \\
未定義 &(x \leq 0, x \geq 1)
\end{cases}
$$
先ほどの形がでてきました。任意の確率密度関数を持つ乱数の生成を行いたい、という場合には上記の手順で実行できます。
理論的背景
なぜ、ある確立密度関数の累積分布関数の逆関数を一様乱数に適用すると、元の確率密度関数に従った乱数が生成されるのか。
一様分布に変換を施す、というのが重要です。ここでは、考えやすくするために、本当に一様な100要素の数列a=[0.00, 0.01, 0.02, 0.03, … 0.99]を考えてみます。小さい順に並んでいます。この要素に対して、何らかの関数を適用した、変換後数列bが特定の確立密度関数に従うものになるようにしたい、ということになります。ここで、”何らかの関数”は、変換前後で要素の大小関係が入れ替わらない（単調増加の）関数とします。
この”何らかの関数”がわからないわけではありますが、数列bの累積分布関数f(x)は、数xが出てくるのは数列bの要素の始めから何％の位置か、を教えてくれる関数と見ることができます。例えば、上の例だと、f(0.25) = 2*0.25-0.25^2 = 0.4375 つまり、約44%の位置（100要素の数列なら約44番目）に0.25が出てくる、ということを教えてくれます。
さて、”何らかの関数”は、一様な数列aを順番を変えずに数列bに変換する関数です。つまり、数列bの要素のうちx%の位置にある数、を教えてくれる関数なのです。つまり、累積分布関数のちょうど逆関数になっているのです。
