確率変数の変数変換は、カルマンフィルタやガウス過程回帰の導出でも登場する、使用頻度の高い手法です。統計検定でも頻出の内容であり、確率変数の和の分布（畳み込み積分）とも密接に関連しています。
数学的にはヤコビアンなどが登場しますが、今回はできるだけイメージができるようにわかりやすく解説していきます。
本記事の内容
- 1次元の変数変換の定理
- ヤコビアンの直感的な意味
- 多次元への拡張
- Python での具体例
前提知識
- 確率密度関数の基本的な性質
- 微分と逆関数の基本
1次元の変数変換
問題設定
確率変数 $X$ の確率密度関数 $f_X(x)$ が既知で、$Y = g(X)$ という変換をしたとき、$Y$ の確率密度関数 $f_Y(y)$ を求めたい、というのが変数変換の問題です。
変数変換の定理
$g$ が単調な関数で逆関数 $g^{-1}$ が存在するとき、
$$ \begin{equation} f_Y(y) = f_X(g^{-1}(y)) \left|\frac{dx}{dy}\right| = f_X(g^{-1}(y)) \left|\frac{d}{dy}g^{-1}(y)\right| \end{equation} $$
直感的な理解
なぜ $\left|\frac{dx}{dy}\right|$ が必要なのでしょうか。
確率密度関数は「確率/区間の幅」の比率を表しています。変数変換により区間の幅が変わるため、確率密度もそれに応じて補正する必要があります。
$x$ の微小区間 $[x, x+dx]$ に含まれる確率は $f_X(x) dx$ です。$y = g(x)$ の変換により、対応する $y$ の区間幅は $|dy| = |g'(x)| |dx|$ になります。確率は保存されるので、
$$ f_Y(y) |dy| = f_X(x) |dx| \implies f_Y(y) = f_X(x) \left|\frac{dx}{dy}\right| $$
具体例1: 線形変換
$X \sim \mathcal{N}(0, 1)$ として、$Y = aX + b$（$a \neq 0$）のとき、
逆変換: $x = (y – b)/a$ なので $\left|\frac{dx}{dy}\right| = \frac{1}{|a|}$
$$ f_Y(y) = f_X\left(\frac{y-b}{a}\right) \cdot \frac{1}{|a|} = \frac{1}{|a|\sqrt{2\pi}} \exp\left(-\frac{(y-b)^2}{2a^2}\right) $$
これは $Y \sim \mathcal{N}(b, a^2)$ であり、期待通りの結果です。
具体例2: 対数正規分布
$X \sim \mathcal{N}(\mu, \sigma^2)$ として、$Y = e^X$ のとき、
逆変換: $x = \ln y$ なので $\left|\frac{dx}{dy}\right| = \frac{1}{y}$（$y > 0$）
$$ \begin{equation} f_Y(y) = \frac{1}{y\sigma\sqrt{2\pi}} \exp\left(-\frac{(\ln y – \mu)^2}{2\sigma^2}\right) \quad (y > 0) \end{equation} $$
これが対数正規分布（Log-Normal Distribution）の確率密度関数です。
多次元の変数変換
$D$ 次元の場合、$\bm{Y} = \bm{g}(\bm{X})$ のとき、
$$ \begin{equation} f_{\bm{Y}}(\bm{y}) = f_{\bm{X}}(\bm{g}^{-1}(\bm{y})) \left|\det\left(\frac{\partial \bm{x}}{\partial \bm{y}}\right)\right| \end{equation} $$
ここで $\frac{\partial \bm{x}}{\partial \bm{y}}$ はヤコビ行列で、その行列式（ヤコビアン）が1次元の $|\frac{dx}{dy}|$ の拡張です。
Box-Muller法
変数変換の有名な応用として、一様分布から正規分布を生成するBox-Muller法があります。
$U_1 \sim U(0,1)$, $U_2 \sim U(0,1)$ として、
$$ \begin{align} X_1 &= \sqrt{-2\ln U_1} \cos(2\pi U_2) \\ X_2 &= \sqrt{-2\ln U_1} \sin(2\pi U_2) \end{align} $$
とすると、$X_1, X_2$ はそれぞれ独立な標準正規分布 $\mathcal{N}(0,1)$ に従います。
Python での実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, lognorm
np.random.seed(42)
n_samples = 100000
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 左上: 線形変換
ax1 = axes[0, 0]
X = np.random.normal(0, 1, n_samples)
a, b = 2.0, 3.0
Y_linear = a * X + b
x_range = np.linspace(-5, 12, 300)
ax1.hist(Y_linear, bins=100, density=True, alpha=0.5, color='blue', label='Y = aX + b (sim)')
ax1.plot(x_range, norm.pdf(x_range, b, abs(a)), 'r-', linewidth=2,
label=f'N({b}, {a**2})')
ax1.set_title(f'Linear: Y = {a}X + {b}', fontsize=12)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
# 右上: 対数正規分布
ax2 = axes[0, 1]
mu_ln, sigma_ln = 1.0, 0.5
X_norm = np.random.normal(mu_ln, sigma_ln, n_samples)
Y_lognorm = np.exp(X_norm)
y_range = np.linspace(0.01, 15, 300)
# 対数正規分布のPDF
pdf_lognorm = (1 / (y_range * sigma_ln * np.sqrt(2*np.pi))) * \
np.exp(-(np.log(y_range) - mu_ln)**2 / (2*sigma_ln**2))
ax2.hist(Y_lognorm, bins=150, density=True, alpha=0.5, color='blue', range=(0, 15), label='Y = exp(X) (sim)')
ax2.plot(y_range, pdf_lognorm, 'r-', linewidth=2, label='Log-Normal (theory)')
ax2.set_title('Log-Normal: Y = exp(X), X ~ N(1, 0.25)', fontsize=12)
ax2.set_xlim(0, 15)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
# 左下: Box-Muller法
ax3 = axes[1, 0]
U1 = np.random.uniform(0, 1, n_samples)
U2 = np.random.uniform(0, 1, n_samples)
X_bm1 = np.sqrt(-2 * np.log(U1)) * np.cos(2 * np.pi * U2)
X_bm2 = np.sqrt(-2 * np.log(U1)) * np.sin(2 * np.pi * U2)
x_range2 = np.linspace(-4, 4, 200)
ax3.hist(X_bm1, bins=100, density=True, alpha=0.4, color='blue', label='X1 (Box-Muller)')
ax3.hist(X_bm2, bins=100, density=True, alpha=0.4, color='green', label='X2 (Box-Muller)')
ax3.plot(x_range2, norm.pdf(x_range2, 0, 1), 'r-', linewidth=2, label='N(0,1) theory')
ax3.set_title('Box-Muller Transform', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
# 右下: ヤコビアンの効果の可視化
ax4 = axes[1, 1]
X_uniform = np.random.uniform(0, 1, n_samples)
Y_sq = X_uniform**2 # Y = X^2
y_sq_range = np.linspace(0.001, 1, 200)
# f_Y(y) = f_X(sqrt(y)) * |dx/dy| = 1 * 1/(2*sqrt(y))
pdf_sq = 1 / (2 * np.sqrt(y_sq_range))
ax4.hist(Y_sq, bins=100, density=True, alpha=0.5, color='blue', label='Y = X^2 (sim)')
ax4.plot(y_sq_range, pdf_sq, 'r-', linewidth=2, label='f_Y(y) = 1/(2*sqrt(y))')
ax4.set_title('Y = X^2 where X ~ U(0,1)', fontsize=12)
ax4.legend(fontsize=10)
ax4.set_xlim(0, 1)
ax4.grid(True, alpha=0.3)
plt.suptitle('Random Variable Transformation', fontsize=14)
plt.tight_layout()
plt.show()
右下のグラフでは、一様分布を二乗する変換によってヤコビアン $\frac{1}{2\sqrt{y}}$ の効果で0付近に密度が集中する様子が確認できます。
まとめ
本記事では、確率変数の変数変換について解説しました。
- 変数変換の定理では、逆変換の微分（ヤコビアン）による補正が必要
- ヤコビアンは変換による区間の幅の変化を補正し、確率を保存する役割を持つ
- 対数正規分布は正規分布の指数変換として導出できる
- Box-Muller法は変数変換を使って一様分布から正規分布を生成する
次のステップとして、以下の記事も参考にしてください。