【統計検定１・準1級】2変量正規分布の条件付き期待値・分散をmatplotlibによる描画で解釈

$X \sim N(\mu_x,\sigma_x^2)$，$Y \sim N(\mu_y,\sigma_y^2)$のとき，
$X=x$が与えられたときの$Y$の条件付き分布$Y|(X=x)$は正規分布となり，その期待値と分散は

\text{条件付き期待値：} E[Y|X=x] = \mu_y + \rho \sigma_y \frac{x-\mu_x}{\sigma_x} = \mu_y + \frac{\sigma_{xy}}{\sigma_x^2}(x-\mu_x)

\text{条件付き分散：} V[Y|X=x] = \sigma_y^2(1-\rho^2)=\sigma_y^2-\frac{\sigma_{xy}}{\sigma_x^2}

となります．

解釈

条件付き期待値の第２項は，「$X$が$\mu_x$からずれていることが，$Y$にどれだけ影響を与えるか」を表しています．

条件付き分散は，分散共分散行列$\Sigma$

\Sigma = 
\begin{pmatrix}
   \sigma_x^2 & \sigma_{xy} \\
   \sigma_{xy} & \sigma_y^2
\end{pmatrix}
=
\begin{pmatrix}
   \sigma_x^2 & \rho \sigma_x \sigma_y \\
   \rho \sigma_x \sigma_y & \sigma_y^2
\end{pmatrix}

を用いて，

\Sigma_{22|1}= \Sigma_{22}-\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}

と表せます．$\Sigma_{11}$は$\Sigma$の$(1,1)$成分を意味しています．
よって，

V[Y|X=x] = \Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{11} = \sigma_x^2-\rho \sigma_x \sigma_y\frac{1}{\sigma_y^2}\rho \sigma_x \sigma_y = \sigma_y^2(1-\rho^2)

条件付き分散$\sigma_y^2(1-\rho^2)$からわかることとして，条件$X=x$が加わることで，$Y$の分散は，条件のない場合の分散$\sigma_y^2$の$(1-\rho^2)$倍になる（小さくなる）のです．

これは，$X$と$Y$の相関が大きいほど，$X$がある値に定まれば$Y$の値のばらつきは小さくなり，$\rho = \pm 1$のとき，$Y$はある値に特定できるため，条件付き分散は０になります．

逆に，$X$と$Y$の相関が小さいほど，条件付き分散はもとの分散に近づき，$\rho=0$のとき，条件付き分散は条件がない場合と同じ$\sigma_y^2$になります．

視覚化

数式だけでは何が起こっているのかわかりにくいので，Pythonを用いて視覚化してみることにします．

2変数$X,Y$がともに正規分布に従い，

\mathbf \mu =
\begin{pmatrix}
   0\\
   0
\end{pmatrix}

\Sigma = 
\begin{pmatrix}
   1 & 0.5 \\
   0.5 & 1
\end{pmatrix}

のとき，2変量正規分布の確率密度関数は以下のようになります．

これを，断面$X=1.0$で切ったときの$(y,z)$平面は次のようになります．

これが，$X=1.0$における条件付き確率密度関数を示しています．
条件付き期待値・分散は，

E[Y|X=1.0] = \mu_y + \frac{\sigma_{xy}}{\sigma_x^2}(x-\mu_x) = 0 + \frac{0.5}{1}(1-0)=0.5

V[Y|X=1.0] =\sigma_y^2-\frac{\sigma_{xy}}{\sigma_x^2} = 1-\frac{0.5}{1}=0.5

と求められます．

元の変数$Y$は，$N(0,1)$に従っていたので，$X=1.0$という条件によって，$N(0.5,0.5)$に従う分布になったことがわかります．

参考までに，上記の描画に使用したコードを載せておきます．（ほとんどchat-gptに丸投げです．）

サンプルコード

pip install numpy matplotlib scipy

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

２変量正規分布を生成し，3Dプロットを作成します．平均meanや分散共分散行列covは，自由に設定してください．

# 平均と分散共分散行列を設定
mean = [0, 0]  # 平均
cov = [[1, 0.5], [0.5, 1]]  # 分散共分散行列

# 2次元グリッドを作成
x, y = np.mgrid[-3:3:.01, -3:3:.01]
pos = np.dstack((x, y))

# 2変量正規分布を定義
rv = multivariate_normal(mean, cov)

# 3Dプロットを作成
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# 確率密度関数を計算
z = rv.pdf(pos)

# 3Dサーフェスプロット
ax.plot_surface(x, y, z, cmap='viridis', edgecolor='none')

# ラベルを追加
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Probability Density')
plt.show()

特定の平面$X=x_{cut}$で分布を切り，その断面を２次元で表示します．

# 切り取るxの値
x_cut = 1.0

# x = x_cutの断面のy方向の確率密度を計算
y_values = np.linspace(-3, 3, 100)
z_cut = rv.pdf(np.column_stack((np.full(100, x_cut), y_values)))

# 断面を2Dプロット
plt.plot(y_values, z_cut, label=f'X = {x_cut}')
plt.xlabel('Y axis')
plt.ylabel('Probability Density')
plt.title(f'Cut at X = {x_cut}')
plt.legend()
plt.grid(True)
plt.show()

この記事がご参考になれば，いいねを押して頂けるととても力になります．