確率変数の変換について - My Life as a Mock Quant

ロジックメモ

ある確率変数 $X: \Omega \rightarrow \mathbb{R}$ とその写像 $Y=g(X): \Omega \rightarrow \mathbb{R}$ で定義される２つの確率変数、特に確率変数 $Y$ の確率分布について考えたい。
$Y$ の確率分布関数 $F_Y(y)$ は
$F_Y(y) = P(Y \le y)= P(g(X) \le y) = P(\left\{ \omega \in \Omega; g(X(\omega)) \le y\right\}) = P(\left\{ \omega \in \Omega; X(\omega) \le g^{-1}(y)\right\}) = P(X \le g^{-1}(y)) = F_X(g^{-1}(y))$
と確率変数 $X$ の確率分布関数を用いてあらわすことができる。ここで、 $g$ は単調増加な関数であると仮定している。単調減少の場合は不等号の向きが逆になる。それ以外の場合は、単調増加・現象とみなせる区間に切って、確率を分けて計算し、最後に足しあげる操作が必要。

したがって、確率変数 $Y$ の確率密度関数 $f_Y$ は確率変数 $X$ の確率密度関数 $f_X$ を用いて
$f_Y(y) = \frac{d}{dy} F_Y(y) = f_X(g^{-1}(y)) \cdot \frac{d}{dy} g^{-1}(y)$
と書くことができる。

Rで確かめる

これを確かめてみよう。ここでは $g(x) = \log x^2$ と定義する。すると上式から
$F_Y(y) = P(\left\{ \omega \in \Omega; \log X(\omega)^2 \le y\right\}) = P(\left\{ \omega \in \Omega; - \exp(y/2) \le X(\omega) \le \exp(y/2) \right\}) = F_X(\exp(y/2)) - F_X(-\exp(y/2))$
とすることができる。ここで関数 $g$ は $x=0$ を境に単調減少から単調増加関数へと変わるので、上のように項が2つ出てくる。上で言及した最後の"足しあげ"に相当するものだ。
従って確率密度関数 $f_Y(y)$ は
$f_Y(y) = f_X(\exp(y/2)) \cdot \frac{1}{2}\exp(y/2) - f_X(-\exp(y/2)) \cdot \frac{1}{2}(-\exp(y/2)) = \frac{1}{2}\exp(y/2) \left\{ f_X(\exp(y/2)) + f_X(-\exp(y/2))\right\}$
となる。更に(本来、何でもいいんだが)確率変数 $X$ が標準正規分布に従うと仮定して上の結果を確かめてみる。

size <- 10^5
#標準席分布に従う乱数の生成
x <- rnorm(size)
#関数gで変換
y <- log(x^2)

この結果をプロットする(確率密度表示)＆上で計算した確率密度関数を上から重ねてみると

hist(y, sqrt(size), freq=FALSE)
ys <- seq(-20, 5, 0.1)
fy <- function(y){exp(y/2)/2*(dnorm(exp(y/2)) + dnorm(-exp(y/2)))}
lines(ys, fy(ys), col=2, lwd=3)

・・・となってほぼぴったり一致する。なので、計算は正しいだろうと。

ggplot2に慣れたいんで、ggplot2でもやってみた。..density..をggplot関数じゃなくてgeom_histogram関数の中に書かなきゃだめなのがミソなんだが、こういうの覚えきれん。。。

library("ggplot2")
ggplot(data.frame(value=y), aes(x=value)) + 
  geom_histogram(aes(y=..density..), binwidth=0.1, fill="light blue") + 
  stat_function(fun=fy, colour="blue", size=2)