Goodな生活

INTP型コンサルタントが好奇心の受け皿を探す

【統計検定準1級】2次元正規分布と条件付き確率分布の求め方【R】

はじめに

この記事では2次元正規分布の密度関数と、2次元のうち1つの確率変数が与えられたときの条件付き確率分布の求め方について扱います。相関係数や回帰分析とも関連するトピックです。

まずは多次元正規分布から考える

2次元正規分布を扱う前に、まずは多次元の正規分布からスタートします。多次元の確率変数は必ずしもそれぞれが独立の確率分布に従う訳ではありません。例えば、確率変数の例としてテストの点数を考えます。ここでは次元の数=教科の数です。生徒にっては文系科目(国語や英語)は低得点だが理系科目(数学や理科)は高得点、というような傾向があると思われます。とすると、算数と理科のテスト点数には何らかの関連がありそうです。このように関連しあう多数の確率変数を仮定する正規分布を、多次元(多変量)正規分布(multidimensional normal distribution)と呼びます。

一般化します。連続確率変数 X_ik次元正規分布の同時密度関数は以下のように表せます。


 {
\begin{eqnarray}
f_{x} = \frac{1}{(2\pi)^{\frac{k}{2}}} \frac{1}{|\Sigma|^{\frac{1}{2}}} exp\left\{-\frac{1}{2}(x -\mu )\top\Sigma^{-1}(x -\mu )\right\}\tag{1}
\end{eqnarray}
}

\muは連続確率変数 X_iの平均、Σ は分散共分散行列(covariance matrix)と呼ばれるk×kの対称行列です。Σ の対角成分には X_iの分散、それ以外は共分散が含まれます。


 {
\begin{eqnarray}
\mu =\begin{pmatrix}
\mu_1 \\
\vdots \\
\mu_k
\end{pmatrix},
\Sigma=\begin{pmatrix}
σ_{11} &\cdots& σ_{1k} \\
\vdots &\ddots& \vdots \\
σ_{k1}  &\cdots & σ_{kk} \\
\end{pmatrix} ,\forall i,j \in {1,...,k} \tag{2}
\end{eqnarray}
}

ここでは\mu = E(X_i), \sigma_{ij} = Cov(X_i,X_j),\sigma_{ii} = V(X_i)です。

2次元正規分布の導出

それでは(1)においてk=2の場合、つまり2次元正規分布の密度関数を導出しましょう。先に指数関数の中の Σ の逆行列を計算します。


 {
\begin{eqnarray}
\Sigma^{-1} &=& \frac{1}{\sigma^{2}_{1}\sigma^{2}_{2}(1-\rho^{2})} 
\begin{pmatrix}
σ^{2}_{2} & -\ \rho σ_{1}σ_{2} \\
 -\ \rho σ_{1}σ_{2}  &  σ^{2}_{1} \\
\end{pmatrix} \\ \tag{3}
\end{eqnarray} 
}

ここでは

{
\begin{eqnarray}
σ_{11} &=& σ^{2}_{1} \\
σ_{11} &=& σ^{2}_{1} \\
\rho &=& \frac{σ_{12} }{σ_{1}σ_{2}} \\
&=& \frac{σ_{21} }{σ_{1}σ_{2}} \tag{4}
\end{eqnarray}
}

ρ相関係数(correlation coefficient)と呼ばれ、標準偏差(分散の平方根をとったもの)と共分散を使って表します。指数関数の中の二次形式を計算すると、確率密度関数は次のような形になります。


 {
\begin{eqnarray}
f_{x} &=&  \frac{1}{2\pi} \frac{1}{\sqrt{1-\rho^2}\sigma_{1}\sigma_{2}} exp\left\{ (x_1 -\mu_1, x_2 -\mu_2) \frac{1}{(1-ρ^2)σ^2_{1}σ^2_{2}} 
\begin{pmatrix}
σ^2_{2} & -σ_{12} \\
\, -σ_{21}  & σ^2_{1} \\
\end{pmatrix}
\begin{pmatrix}
x_1 -\mu_1 \\
x_2 -\mu_2 \\
\end{pmatrix}
\right\} \\
&=& \frac{1}{2\pi} \frac{1}{\sqrt{1-\rho^2}\sigma_{1}\sigma_{2}} 
exp\left\{-\frac{1}{2(1-\rho^2)}\left[\left(\frac{x_1-\mu_1}{σ_1}\right)^2  - 2\rho \frac{x_1-\mu_1}{σ_1}\frac{x_2-\mu_2}{σ_2}  + \left(\frac{x_2-\mu_2}{σ_2}\right)^2\right]\right\} \tag{5}\\
\end{eqnarray}
}


指数関数の内側はX_1,X_2を標準化した変数の楕円の方程式になっており、相関係数ρの大きさによって楕円の形が変化します。ρが大きくなるほど、楕円の形がよりシャープに(線形に近づく)なります。具体的な値をρに代入して、楕円の振る舞いがどのように変化するかを確認しましょう。

グラフでみる2次元正規分布

相関係数ρ=0.5, 0.9のそれぞれの場合の、2次元正規分布から得られたデータ(実現値)の散布図です。ここでは単純化のためσ_1=σ_2=1と仮定します。

ρ=0.5とρ=0.9の散布図の比較

相関係数の値が大きくなる=線形な関係性が強くなることが分かります。
同様に相関係数ρ=0.5, 0.9のそれぞれの場合の、2次元正規分布の確率密度関数を示します。

ρ=0.5の密度関数
ρ=0.9の密度関数

相関係数ρの値が大きくなるにつれて、2変数の関係が強くなり、密度関数の形がよりシャープになる様子が分かります。

2次元正規分布を使って「無相関なら独立」を示す

2つの確率変数に相関関係が存在しないとき、つまりρ=0のとき、2次元正規分布の密度関数はどのような形になるのでしょうか。

(5)にρ=0を代入すると、


 {
\begin{eqnarray}
f_{x} &=& \frac{1}{2\pi} \frac{1}{\sigma_{1}\sigma_{2}} 
exp\left\{-\frac{1}{2}
\left[
\left(\frac{x_1-\mu_1}{σ_1}\right)^2 + \left(\frac{x_2-\mu_2}{σ_2}\right)^2
\right]
\right\} \\

&=& \frac{1}{\sqrt{2\pi\sigma^{2}_{1}}}  
exp\left\{-\frac{1}{2}
\left(\frac{x_1-\mu_1}{σ_1}\right)^2 
\right\}

\frac{1}{\sqrt{2\pi\sigma^{2}_{2}}}
exp\left\{-\frac{1}{2}
 \left(\frac{x_2-\mu_2}{σ_2}\right)^2
\right\}\tag{6} 
 
\end{eqnarray}
}

(6)では、2次元正規分布の密度関数が2つの密度関数の積で表されています。これは2つの確率変数が無相関ならば独立であるという2次元正規分布に従う確率変数の重要な性質です。

条件付き確率分布の求め方 (パターン1)

最後に、2次元正規分布を使った条件付き確率分布(密度関数)の求め方を紹介します。条件付き確率分布とは、2つの確率変数のうち1つの確率変数が与えられたときの、もう片方の確率変数の条件付き分布です。ここでは2次元確率変数(X_1,X_2)のうちX_2 = x_2が与えられたときのX_1の分布を考えます。求め方は2パターンあります。まず密度関数を変形する方法です。

(5)の指数関数の中身を、x_1でまとめた項と、その残りに分けます。


 {
\begin{eqnarray}
f_{x}
&=& \frac{1}{2\pi} \frac{1}{\sqrt{1-\rho^2}\sigma_{1}\sigma_{2}} 
exp\left\{-\frac{1}{2(1-\rho^2)}\left[\left(\frac{x_1-\mu_1}{σ_1}\right)^2  - 2\rho \frac{x_1-\mu_1}{σ_1}\frac{x_2-\mu_2}{σ_2}  + \left(\frac{x_2-\mu_2}{σ_2}\right)^2\right]\right\} \\

&=& \frac{1}{2\pi} \frac{1}{\sqrt{1-\rho^2}\sigma_{1}\sigma_{2}}
exp\left\{-\frac{1}{2(1-\rho^2)}\left[
\frac{1}{σ^2_1} \left(x_1 - \mu_1 - \rho \frac{σ_1}{σ_2} (x_2 -\mu_2)\right)^2  + \frac{1-\rho^2}{σ^2_2} \left(x_2 - \mu_2\right)^2
\right]
\right\} \\

&=& \frac{1}{2\pi} \frac{1}{\sqrt{1-\rho^2}\sigma_{1}\sigma_{2}}
exp\left\{-\frac{1}{2(1-\rho^2)}\frac{1}{σ^2_1} \left(x_1 - \mu_1 - \rho \frac{σ_1}{σ_2} (x_2 -\mu_2)\right)^2  + \frac{1}{σ^2_2} \left(x_2 - \mu_2\right)^2
\right\}\tag{7}\\
\end{eqnarray}
}

(7)より、X_1の条件付き期待値E[X_1|X_2=x_2]\mu_1 + \rho \frac{σ_1}{σ_2}(x_2 -\mu_2)です。
これはX_1X_2に回帰させる単回帰の式を使っても表すことができます。

{
\begin{eqnarray}
E[X_1|X_2=x_2] &=& \mu_1 + \rho \frac{σ_1}{σ_2}(x_2 -\mu_2) \\
                           &=& \mu_1 -  \rho \frac{σ_1}{σ_2}\mu_2 + \rho \frac{σ_1}{σ_2} x_2 \\
                           &=& \bar{X_1} - \hat{β_2} \bar{X_2} +  \hat{β_2} x_2 \\
                           &=& \hat{β_1} +  \hat{β_2} x_2 \tag{8}
\end{eqnarray}
}


(8)の\hat{β_1},\hat{β_2}は以下の単回帰の最小二乗推定量です。

{
\begin{eqnarray}
X_1 &=& β_1 + β_2 X_2 + ε \tag{9} 
\end{eqnarray}
}


(8)が意味するのは、X_2=x_2が与えられたときのX_1の条件付き期待値は、単回帰によって導出された予測値と等しくなるという性質です。
同じく(7)より、条件付き分散は、σ^2_1(1-ρ^2)です。X_2が所与(非確率変数)であるため、X_2のばらつきであるσ_2は登場しません。

条件付き確率分布の求め方 (パターン2)

2つ目のパターンは、確率変数(X_1,X_2)をそれぞれ標準化した確率変数を用いて、条件付きの平均と分散を求める方法です。

(X_1,X_2)をそれぞれ標準化すると、

{
\begin{eqnarray}
X^{*}_1 = \frac{X_1 -μ_1}{σ_1} , \, X^{*}_2 = \frac{X_2 - μ_2}{σ_2} \tag{10}
\end{eqnarray}
}

 {
\begin{eqnarray}
\begin{pmatrix}
X^{*}_1 \\
X^{*}_2 
\end{pmatrix}
\rightarrow N
\left(
\begin{pmatrix}
0 \\
0
\end{pmatrix},
\begin{pmatrix}
1 & ρ \\
ρ & 1
\end{pmatrix}
\right) \tag{11}
\end{eqnarray}
}


となるため、Z = X^{*}_1 - ρX^{*}_2とおけば、Zの平均はE[Z]=0、分散は


{
\begin{eqnarray}
V[Z] &=& V(X^{*}_1) - 2ρ Cov(X^{*}_1,X^{*}_2)+ρ^2 V(X^{*}_2) \\
           &=& 1 - ρ^2 \tag{12}
\end{eqnarray}
}

です。ここで(10)より、

{
\begin{eqnarray}
X_1 = μ_1 + X^{*}_1 σ_1 \tag{13}
\end{eqnarray}
}


です。(13)の両辺をX_2=x_2で条件付けます。X_1,X^{*}_1はどちらも確率変数であるため、条件付き期待値が両辺に登場します。X_1の平均と分散は、


{
\begin{eqnarray}
E[X_1|X_2=x_2] &=& μ_1 + σ_1 E[X^{*}_1|X_2=x_2] \\
         &=& μ_1 + σ_1 E[Z + ρ X^{*}_2|X_2=x_2] \\
                           &=& μ_1 + σ_1 E[Z + ρ\frac{x_2 -μ_2}{σ_1}] \\
                           &=& μ_1 + ρ σ_1 \frac{x_2-μ_2}{σ_2} \tag{14}
\end{eqnarray}
}

{
\begin{eqnarray}
V[X_1|X_2 = x_2] &=& σ^2_{1} V[X^{*}_1 | X_2 = x_2]  \\
                              &=& σ^2_{1} V[Z + ρ x_2] \\
                              &=& σ^2_{1} (1-ρ^2) \tag{15}
\end{eqnarray}
}

と求めることができました。

参考文献


確率密度関数を書くためのRのソースコードは三中信宏教授の作成されたものを参考にさせていただきました。
http://cse.naro.affrc.go.jp/minaka/R/R-binormal.html

読んでいただいてありがとうございました。