Zenn
💨

二変量正規分布においてaX₁+bX₂が従う分布

2022/07/13に公開
Python
統計学
正規分布
tech

求めるもの

\begin{pmatrix} X_1 \\ X_2 \\ \end{pmatrix} \sim \mathcal{N}(\mathbf{\mu}, \mathbf{\Sigma}), \mathbf{\mu} = \begin{pmatrix} \mu_1 \\ \mu_2 \\ \end{pmatrix}, \mathbf{\Sigma} = \begin{pmatrix} \sigma_{11} & \sigma_{12} \\ \sigma_{21} & \sigma_{22} \\ \end{pmatrix}

上記のように、X_1, X_2が二変量の正規分布に従う時a, bを実数として、a X_1 + b X_2が従う分布を考えます。

多変量正規分布の線形変換に関する定理

\mathbf{X} \sim \mathcal{N}(\mathbf{\mu}, \mathbf{\Sigma})のとき、正則行列Cに対して、

\begin{aligned} C\mathbf{X} \sim \mathcal{N}(C\mathbf{\mu}, C\mathbf{\Sigma}C^T) \end{aligned}

導出

\begin{pmatrix} Y_1 \\ Y_2 \\ \end{pmatrix} = C\begin{pmatrix} X_1 \\ X_2 \\ \end{pmatrix}, C = \begin{pmatrix} 1 & -\frac{\sigma_{12}}{\sigma_{22}} \\ 0 & 1 \\ \end{pmatrix}

という変換を考えます。ここで、Cは常に正則行列であるので、多変量正規分布の線形変換に関する定理を用います。

\begin{aligned} \begin{pmatrix} Y_1 \\ Y_2 \\ \end{pmatrix} &\sim \mathcal{N}(C\mathbf{\mu}, C\mathbf{\Sigma}C^T) \\ where \quad C\mathbf{\mu} &= \begin{pmatrix} \mu_1 - \frac{\sigma_{12}}{\sigma_{22}}\mu_2 \\ \mu_2 \\ \end{pmatrix} \\ C\mathbf{\Sigma}C^T &= \begin{pmatrix} \sigma_{11} - \frac{\sigma_{12}^2}{\sigma_{22}} & 0 \\ 0 & \sigma_{22} \\ \end{pmatrix} \end{aligned}

今、分散共分散行列が対角行列となっているので、

\begin{aligned} Y_1 &\sim \mathcal{N}(\mu_1 - \frac{\sigma_{12}}{\sigma_{22}}\mu_2, \sigma_{11} - \frac{\sigma_{12}^2}{\sigma_{22}}) \\ Y_2 &\sim \mathcal{N}(\mu_2, \sigma_{22}) \end{aligned}

のように表すことができます。また、Y_1, Y_2は互いに独立です。以上を用いて、

\begin{aligned} aX_1 + bX_2 &= aY_1 + (b + a\frac{\sigma_{12}}{\sigma_{22}})Y_2 \\ & \sim \mathcal{N}(\mathbf{\mu}^*, \mathbf{\sigma^*}) \\ where \quad \mathbf{\mu}^* &= a(\mu_1 - \frac{\sigma_{12}}{\sigma_{22}}\mu_2) + (b + a\frac{\sigma_{12}}{\sigma_{22}})\mu_2 \\ \mathbf{\sigma^*} &= a^2(\sigma_{11} - \frac{\sigma_{12}^2}{\sigma_{22}}) + (b + a\frac{\sigma_{12}}{\sigma_{22}})^2\sigma_{22} \end{aligned}

となります。

コードで確認

結果が正しいかどうかをPythonを用いて確認しました。青い線は、\mathcal{N}(\mathbf{\mu}^*, \mathbf{\sigma^*})の確率密度関数です。オレンジの線は、\mathbf{X} \sim \mathcal{N}(\mathbf{\mu}, \mathbf{\Sigma})から乱数により10^6個のデータをサンプリングし、aX_1 + bX_2の値を求めたものです。

pythonコード
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import multivariate_normal, norm

np.random.seed(28)

fig = plt.figure(figsize=(12, 20))

for i in range(10):
    a, b = np.random.standard_cauchy(2)
    mu = np.random.standard_cauchy(2)
    sigma_11, sigma_22 = abs(np.random.standard_cauchy(2))
    sigma_12 = np.sqrt(sigma_11) * np.sqrt(sigma_22) * (np.random.rand() - 0.5) * 2
    sigma = np.array([[sigma_11, sigma_12],
                      [sigma_12, sigma_22]])

    data = multivariate_normal(mu, sigma).rvs(size=int(1e6))

    x1, x2 = data[:, 0], data[:, 1]
    ax1_bx2 = a * data[:, 0] + b * data[:, 1]

    mu_star = a*(mu[0] - sigma_12/sigma_22*mu[1]) + (b + a*sigma_12/sigma_22)*mu[1]
    sigma_star = a**2*(sigma_11 - sigma_12**2/sigma_22) + (b + a*sigma_12/sigma_22)**2*sigma_22

    xs = np.linspace(
        mu_star - 3*np.sqrt(sigma_star),
        mu_star + 3*np.sqrt(sigma_star),
        num=1000
    )
    dens = norm.pdf(
        x=xs, loc=mu_star, scale=np.sqrt(sigma_star)
    )
    data_bar = [0 for _ in range(len(xs))]
    for x in ax1_bx2:
        ind = (x - mu_star + 3*np.sqrt(sigma_star) + 6*np.sqrt(sigma_star)/len(xs)/2) * len(xs) / (6 * np.sqrt(sigma_star))
        ind = int(ind)
        if (0 <= ind < len(xs)):
            data_bar[ind] += 1
    
    ax = fig.add_subplot(5, 2, i+1)
    ax.plot(xs, dens, label="N({:.2f}, {:.2f})".format(mu_star, sigma_star))
    ax.plot(xs, np.array(data_bar) / sum(data_bar) * (len(xs) / (6 * np.sqrt(sigma_star))), alpha=0.5, label="sample data")
    ax.legend(loc="upper right")

    print("mu_star: {:.2f}, sigma_star: {:.2f}".format(mu_star, sigma_star), end=" ")
    print("ax1_bx2の平均: {:.2f}, ax1_bx2の分散: {:.2f}".format(ax1_bx2.mean(), ax1_bx2.var()))

plt.show()

mu_star: -0.29, sigma_star: 0.10 ax1_bx2の平均: -0.29, ax1_bx2の分散: 0.10
mu_star: 0.25, sigma_star: 1.07 ax1_bx2の平均: 0.25, ax1_bx2の分散: 1.07
mu_star: 2.31, sigma_star: 2.05 ax1_bx2の平均: 2.31, ax1_bx2の分散: 2.05
mu_star: 0.00, sigma_star: 28.00 ax1_bx2の平均: 0.00, ax1_bx2の分散: 27.98
mu_star: -57.62, sigma_star: 23768.73 ax1_bx2の平均: -57.47, ax1_bx2の分散: 23759.03
mu_star: 3.29, sigma_star: 0.58 ax1_bx2の平均: 3.29, ax1_bx2の分散: 0.58
mu_star: -4.42, sigma_star: 23.98 ax1_bx2の平均: -4.42, ax1_bx2の分散: 23.98
mu_star: -4.38, sigma_star: 5.91 ax1_bx2の平均: -4.38, ax1_bx2の分散: 5.91
mu_star: -3.02, sigma_star: 114.45 ax1_bx2の平均: -3.03, ax1_bx2の分散: 114.38
mu_star: -17.03, sigma_star: 211.29 ax1_bx2の平均: -17.02, ax1_bx2の分散: 211.23

density_and_sample

Discussion