2次元確率分布の変数変換（計算してPythonで確認）

この記事では2次元確率分布の変数変換を紹介します。
ヤコビアンを使うやつです。

はじめに公式を確認してから例題で実際に計算し、最後にPythonで確認してみます。
自分が学習した備忘録ですが、2次元確率分布の変数変換の勉強をしているけどよく分からない、イメージが湧かない、という方に助けになれば嬉しいです。

変数変換とは？

2次元確率分布の変数変換とは、確率変数$\ (X, Y)$の同時確率（密度）関数を$\ f_{X,Y}(x, y)$ とするとき、$\ S = g_1(X, Y), T = g_2(X, Y)$ という変換を施すことです。

そして大抵は 変換後の$(S, T)$の同時確率（密度）関数$\ f_{S,T}(s, t)$ がどんな分布になるのかを知りたい！ となります。
この記事ではそれを知るための計算方法を紹介します。

（1次元の場合ですが）変数変換の具体例として下記のようなものがあります。

• データの標準化
  正規分布${\it N}(\mu, \sigma^2)$に従う確率変数$X$に

  $$Z=(X - \mu)/\sigma$$

  なる変数変換を施すと、確率変数$Z$は標準正規分布${\it N}(0, 1)$に従います。
  データの標準化と言われるやつですね。

• 時間の単位の変換
  もっと単純な例ですが、分単位で取られた何かの計測時間データに60を掛けて秒単位のデータに変えるのも変数変換です。
  この場合は、分単位データを$X$とすると、$Y = X \times 60$ の変数変換となります。
  もちろんこの変換では分布の形状は変わりません。

このような変換を2次元で行い、変換後の分布を知りたいというのがモチベーションです。

変数変換の公式

2次元確率分布の変数変換の公式は下記です。

確率変数$\ (X, Y)$の同時確率（密度）関数を$\ f_{X,Y}(x, y)$ とし、
$X = h_1(S, T), Y = h_2(S, T)$と表されるとき、$(S, T)$の同時確率（密度）関数$f_{S,T}(s, t)$は

$$f_{S,T}(s,t) = f_{X,Y}(h_1(s, t), h_2(s, t)) |J(s,t)|$$

となる。

式中の$J(s,t)$は ヤコビアン といい、

$$J(s,t) = det \left( \begin {array} {cc} (\partial/\partial s) h_1(s,t) & (\partial/\partial t) h_1(s,t) \\ (\partial/\partial s) h_2(s,t) & (\partial/\partial t) h_2(s,t) \end {array} \right)$$

で求めます。ただし、行列${\bold A}$に対して、$det({\bold A})$は行列式を表し

$$det \left( \begin {array} {cc} a & b \\ c & d \end {array} \right) = ad - bc$$

となります。

この記事では上記公式やヤコビアンについて詳しく説明しませんが、以降で例題を解いてみます。

変数変換の例

確率変数$X, Y$が独立に分布していてそれぞれガンマ分布に従い、

$$X \sim {\it Ga(a, 1)} \\ Y \sim {\it Ga(b, 1)}$$

であるとき、（ただし、$a > 0, b > 0$）

$$S = X + Y \\ T = X / (X + Y)$$

という変数変換をすると、$(S, T)$の同時確率（密度）関数$f_{S,T}(s, t)$はどうなるか実際に計算して求めてみましょう。

先に結論を書いてしまうと、$S, T$は独立に分布し、それぞれガンマ分布とベータ分布に従い

$$S \sim {\it Ga}(a+b,1) \\ T \sim {\it Beta}(a, b)$$

となります。
$S,T$は独立なので$f_{S,T}(s,t) = f_S(s)f_T(t)$となります。
実際に計算してこうなることを確認してみましょう。

---

ガンマ分布とベータ分布

まずこの例題に登場するガンマ分布とベータ分布について確認します。

ガンマ分布

ガンマ分布の確率密度関数は下記です。

$${\it Ga}(\alpha, \beta) = f(x | \alpha, \beta) = \frac{1}{\Gamma(\alpha) \beta^\alpha} x^{\alpha-1} e^{-x/\beta} \ \ , \ \ x > 0$$

ただし、$\alpha > 0, \beta > 0$ となります。

式中の$\ \Gamma(\alpha)$ はガンマ関数

$$\Gamma(\alpha) = \int^\infty_0 y^{\alpha -1} e^{-y} dy$$

です。


Pythonでどんな分布か確認してみましょう。

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

plt.style.use("seaborn")


def gamma_func(a):
    """ガンマ関数"""
    y, abserr = integrate.quad(lambda x: x**(a - 1) * np.exp(-x), 0, np.inf)
    return y

def gamma_pdf(x, a, b):
    """ガンマ分布の確率密度関数"""
    return 1 / (gamma_func(a) * b**a) * x**(a-1) * np.exp(-x / b)


X = np.linspace(0, 10, 1000)

fig, ax = plt.subplots(figsize=(10, 8))
for a, b in zip([1, 2, 4], [1, 1, 1]):
    ax.plot(X, gamma_pdf(X, a, b), label=f"$Ga({a}, {b})$")
    
ax.set(xlabel="X")
ax.legend(fontsize=15)
plt.show()

ちなみに指数分布$\ {\it Ex}(\lambda)$はガンマ分布の特別な場合であり、${\it Ga}(1, 1/\lambda)$と表すことができます。

ベータ分布

ベータ分布の確率密度関数は下記です。

$${\it Beta}(a, b) = f(x | a,b) = \frac{1}{B(a,b)} x^{a-1} (1-x)^{b-1} \ \ , \ \ 0 < x < 1$$

ただし、$a>0, b>0$ となります。

式中の$\ B(a,b)$ はベータ関数

$$B(a,b) = \int^1_0 x^{a-1} (1-x)^{b-1} dx$$

です。


ベータ分布もPythonで確認してみましょう。

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

plt.style.use("seaborn")


def beta_func(a, b):
    """ベータ関数"""
    y, abserr = integrate.quad(lambda x: x**(a-1) * (1 - x)**(b-1), 0, 1)
    return y


def beta_pdf(x, a, b):
    """ベータ分布の確率密度関数"""
    return (1 / beta_func(a, b)) * x**(a-1) * (1 - x)**(b-1)


X = np.linspace(0, 1, 1000)

fig, ax = plt.subplots(figsize=(10, 8))

for a, b in zip([0.1, 1, 3], [0.5, 1, 2]):
    ax.plot(X, beta_pdf(X, a, b), label=f"$Beta({a}, {b})$")
    
ax.set(xlabel="X", ylim=(-0.1, 4))
ax.legend(fontsize=15)
plt.show()

なかなか面白い形をしていますね。
ベータ分布は$a,b$の取り方で色々な形にすることができます。

ガンマ関数とベータ関数の関係

ガンマ関数とベータ関数の間には下記の関係があり、今回の例題を解くのに必要となってきます。

$$B(a,b) = \frac{\Gamma(a) \Gamma(b)}{\Gamma(a+b)}$$

---

計算してみる

ガンマ分布とベータ分布の確認でやや脱線しましたが、例題を実際に計算して解いてみましょう。

まずヤコビアンを求めます。

$$S = X + Y \\ T = X / (X + Y)$$

なので

$$x = h_1(s,t) = st \\ y = h_2(s,t) = s(1 - t)$$

となります。
よってヤコビアンは

$$J(s,t) = det \left( \begin {array} {cc} t & s \\ 1-t & -s \end {array} \right) = -st - (1-t)s = -s$$

となります。

これを変数変換の公式に入れると

$$\begin{aligned} f_{S,T}(s,t) &= f_{X,Y}(h_1(s,t), h_2(s,t))\ |J(s,t)| \\ &= f_{X,Y}(st,s(1-t)) s \\ &= f_X(st)f_Y(s(1-t)) s \\ &= \frac{1}{\Gamma(a)} \{st\}^{a-1} e^{-st} \frac{1}{\Gamma(b)} \{s(1-t)\}^{b-1} e^{-s(1-t)} s \\ &= \frac{1}{\Gamma(a) \Gamma(b)} \{st\}^{a-1} \{s(1-t)\}^{b-1} e^{-s} s \\ &= \frac{1}{\Gamma(a) \Gamma(b)} s^{a+b-1} t^{a-1} (1-t)^{b-1} e^{-s} \\ &= \frac{1}{\Gamma(a+b)} s^{a+b-1} e^{-s} \times \frac{\Gamma(a+b)}{\Gamma(a) \Gamma(b)} t^{a-1} (1-t)^{b-1} \\ \end{aligned}$$

と書けます。ただし$s>0,\ 0<t<1$ です。

info:
式の3行目は$X,Y$が独立なので$f_{X,Y}(x,y) = f_X(x)f_Y(y)$と変形しています。
7行目は6行目の式に$\ \Gamma(a+b) / \Gamma(a+b)$ を掛けて並べ替えています。

さて、式の最後の行をよく見ると「$\times$」の左側は$\ S \sim {\it Ga}(a+b,1)$ を示しています。
また「$\times$」の右側は$\ \Gamma(a+b) / \Gamma(a) \Gamma(b) = 1 / B(a, b)$ に注意すると$\ T \sim {\it Beta}(a,b)$ を示していることが分かります。

したがって、求めたかった$S,T$の確率密度関数$\ f_{S,T}(s,t)$ はガンマ分布とベータ分布と積

$$\begin{aligned} f_{S,T}(s,t) &= \frac{1}{\Gamma(a+b)} s^{a+b-1} e^{-s} \times \frac{1}{B(a,b)} t^{a-1} (1-t)^{b-1} \\ &= {\it Ga}(a+b,1) {\it Beta}(a,b) \end{aligned}$$

と書けます。これはガンマ分布とベータ分布の単純な積なので$S,T$が独立に分布していることも示しています。

まとめると、

確率変数$\ X,Y$ が独立に分布していてそれぞれガンマ分布に従い、

$$X \sim {\it Ga(a, 1)}\\ Y \sim {\it Ga(b, 1)}$$

であるとき、（ただし、$a > 0, b > 0$）

$$S = X + Y \\ T = X / (X + Y)$$

という変数変換をすると、$S, T$は独立に分布し、それぞれガンマ分布とベータ分布

$$S \sim {\it Ga}(a+b,1) \\ T \sim {\it Beta}(a, b)$$

に従う。

ということです！

---

Pythonで確認してみる

最後に上記の変数変換が本当にそうなるのか確認してみましょう。
今回は$\ a=3, b=2$ として確率変数$\ X,Y$ を生成し、変数変換します。

まず、$\ X,Y$ を生成して分布を確認してみましょう。

# 変数変換前の確率変数X,Yの生成
a = 3
b = 2

X = np.linspace(0, 20, 1000)
Y = np.linspace(0, 20, 1000)

# 乱数生成
x = np.random.gamma(a, 1, 100_000)
y = np.random.gamma(b, 1, 100_000)


fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6), dpi=80)

# 上で作成した関数 gamma_pdf()を使っているのでエラーに注意.
ax1.hist(x, bins=50, density=True)
ax1.plot(X, gamma_pdf(X, a, 1), label=f"$Ga({a}, 1)$", c="r")
ax2.hist(y, bins=50, density=True)
ax2.plot(Y, gamma_pdf(Y, b, 1), label=f"$Ga({b}, 1)$", c="r")

ax1.set_xlabel("X", fontsize=13)
ax2.set_xlabel("Y", fontsize=13)
ax1.legend(fontsize=15)
ax2.legend(fontsize=15);

$X,Y$ が$\ {\it Ga}(3,1), {\it Ga}(2,1)$ に従っていることが確認できました。
続いて変数変換

$$S = X + Y \\ T = X / (X + Y)$$

を行います。
先ほどの計算が正しければ$\ S,T$はそれぞれ$\ {\it Ga}(5,1), {\it Beta}(3,2)$ に従っているはずです。

# 変数変換後の乱数が本当に上記のようになっているか確認

S = np.linspace(0, 20, 100_000)
T = np.linspace(0, 1, 100_000)

# 変数変換
s = x + y
t = x / (x + y)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 12), dpi=80)

# 上で作成した関数 gamma_pdf(), beta_pdf()を使っているのでエラーに注意.
ax1.hist(s, bins=50, density=True)
ax1.plot(S, gamma_pdf(S, a + b, 1), label=f"$Ga({a + b}, 1)$", color="r")
ax2.hist(t, bins=50, density=True)
ax2.plot(T, beta_pdf(T, a, b), label=f"$Beta({a}, {b})$", color="r")

ax1.set_title("S", fontsize=15)
ax2.set_title("T", fontsize=15)
ax1.legend(fontsize=15)
ax2.legend(fontsize=15);

${\it Ga}(5,1), {\it Beta}(3,2)$に従っていそうです！

---

というわけでここまで2次元確率分布の変数変換を見てきました。
実際に手を動かして計算してみるとそれほど大変ではないかなと思いますが、ただ計算しているだけだとイメージが湧きづらいので今回はPythonで確認まで行ってみました。

こういった記事は初めて書くので、分かりづらい点や間違い等あるかもしれませんが、何かありましたらコメント等でご教授いただけると幸いです！