確率密度関数 \(f\) の累積分布関数の逆関数が簡単に求まるのであれば、確率変数 \(U \sim U(0, 1)\) を使って分布\(f\)に従う変数を作ることができます。

実際、 \(F(x) = P(X \leq x) = \int_{-\infty}^{x}f(x)\mathrm{d}x\) の逆関数 \(F^{-1}(x)\) が存在した時、 \(X = F^{-1}(U)\) とおくと、 \(F\) の単調性より \(P(X \leq x) = P(F^{-1}(U) \leq F^{-1}(u)) = P(U \leq u)\) となり、これを計算すると、

\[P(U \leq u) = \int_{0}^{u}\mathrm{d}u = u\] \[= F(x) = \int_{-\infty}^{x}f(x)\mathrm{d}x\]

となり、\(X\)の分布が\(f\)であることがわかります。

例として、指数分布\(\mathrm{Ex}(\lambda)\)に従う乱数を生成してみます。

指数分布は以下の確率密度関数で表される分布です。

\[f_{\mathrm{Ex}(\lambda)}(x) = \lambda e^{-\lambda x}\]

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

\[F(x) = P(x \leq X) = \int_{0}^{x}f_{\mathrm{Ex}(\lambda)}(x)\mathrm{d}x\] \[= [-e^{-\lambda x}]_{0}^{x} = 1 - e^{-\lambda x}\]

この逆関数を求めてみると、以下のようになります。

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

というわけで、\(U \sim U(0, 1)\)のとき、\(X = F^{-1}(U) = -\frac{\log(1 - u)}{\lambda}\)が指数分布となることがわかりました。

これを実際に実装してみます。一様乱数は生成できるものとして、そこから逆関数法による指数分布を作り、Rの組み込み関数rexp()の生成する乱数と比較してみました。

library(ggplot2)
library(tidyr)
library(dplyr)

N <- 10000
df <- data.frame(unif = runif(N, min = 0, max = 1)) # 一様乱数

lambda <- 1
df$exp1 <- rexp(N, rate = lambda) # 指数分布(標準)
df$exp2 <- -log(1 - df$unif) / lambda # 指数分布(自作)

df.plot <- df %>% select(one_of(c("exp1", "exp2"))) %>% gather()
ggplot(df.plot, aes(x = value, fill = key))
  + geom_histogram(bins = 10, boundary = 1, position = "identity", alpha = 0.5)

/img/post/2017-07-22-rexp-hist.png

exp1rexp()の生成する乱数で、exp2が独自のものです。

両者がほとんど一致していることがわかります。