逆変換法

無作為な観測データ$X$を生成するために、最初に一様分布の確率変数$u$を作成し、$F^{-1}_{X}(u)$の値に逆変換する。$F^{-1}_{X}(u)$を計算できるのであれば、この方法は利用可能である。つまり、下記の手順で実行すれば変換できる。

逆関数$F^{-1}_{X}(u)$を作成
$F^{-1}_{X}(u)$を計算する
必要な確率変数について、一様分布から確率変数$u$を生成
$x = F^{-1}_{X}(u)$を与える

ここでは、確率密度関数$f_{X}(x)=3x^{2}(0 < x <1)$で確率分布から無作為標本をシュミレーションするために逆変換を行う。確率密度関数$f_{X}(x)=3x^{2}$を積分すると、累積分布関数$ F_{X}(x)=x^3 , (0 < x < 1)$およびその逆関数の$F^{-1}_{X}(u)=u^{1/3}$が得られる。ベクトル$u$として一様乱数を生成する。

言い換えると、$f_{X}(x)=3x^{2}$の確率密度関数に従うような乱数が欲しい場合、その確率密度関数の累積分布関数の逆関数を作り、それに一様分布に従う乱数を入力すれば、$f_{X}(x)=3x^{2}$の確率密度関数に従う乱数が生成できる。

逆変換法によって生成した無作為標本の確率密度ヒストグラムと理論的密度曲線のプロットを表示する。

n <- 1000
u <- runif(n)
# 累積分布関数F(X)の逆関数F^-1(u)
x <- u^(1/3)
hist(x, prob = TRUE, main = bquote('f(x)=3x^2'))
y <- seq(0, 1, 0.01)
lines(y, 3*y^2, lwd = 2, xlab = "", ylab = "", col = 'tomato')

期待値

この話の続きとして、$X$の期待値を考える。つまり、$F_{X}^{-1}(U)$の平均値を考える。$U$は0と1の間を一様にとる。そこで、$X$が確率変数であるとき、分布関数の逆関数を用いて、$X$の期待値$E[X]$を定義する。

\[ E[X] = \int_{0}^{1} F_{X}^{-1}(u) du \]

分布関数の逆関数によって積分の結果は決まるので、「同じ分布に従う確率変数の期待値は等しい」といえる。

数値計算では積分が発散せずに何らかの値が出る場合もあるので、関数$g$の定義域が$X$の取りうる値の範囲を含む場合、確率変数$g(X)$の期待値$E[g(X)]$は下記のように定義される。

\[ E[g(X)] = \int_{0}^{1} g(F_{X}^{-1}(u)) du \]

つまり下記のように、標準正規分布に従い、自然対数をとられているものの値を期待値を知りたいとする。

set.seed(1989)
x <- exp(rnorm(n = 10000, mean = 0, sd = 1))
# これがわからないとする
mean(x)

## [1] 1.642276

確率変数の分布関数の逆関数と一様乱数を使って0から1まで積分すれば、期待値を計算できる。標準正規分布の累積分布関数の逆関数はqnorm()で計算できる。つまり、標準正規分布の分布関数の累積分布関数の自然対数をとった関数を0から1まで積分した結果が期待値となる。

exp_normal_func <- function(p) exp(qnorm(p, mean = 0, sd = 1))
integrate(f = exp_normal_func, lower = 0, upper = 1)

## 1.648723 with absolute error < 8.4e-05

まだまだ統計学の世界は知らないことばかり・・・。