2変量正規分布の等高線（確率楕円）の描き方 - ジョンとヨーコのイマジン日記

計算ノート

2変量正規分布の密度関数 $f(\boldsymbol{x})$ は次式で表される：

$\displaystyle f(\boldsymbol{x}) = \frac{1}{2\pi \det( \Sigma ) }\exp\left(-{1 \over 2} (\boldsymbol{x}-\boldsymbol\mu)^{\top}\Sigma^{-1} ({\boldsymbol x}-\boldsymbol \mu)\right)$

ここで $\boldsymbol{x}$ は2次元のベクトル.

分散共分散行列の逆行列を $\Lambda = \Sigma^{-1}$ と置くと,

$\displaystyle f(\boldsymbol{x}) = \frac{\det( \Lambda )}{2\pi}\exp\left(-{1 \over 2} (\boldsymbol{x}-\boldsymbol\mu)^{\top} \Lambda ({\boldsymbol x}-\boldsymbol \mu)\right)$

と書ける.

密度関数の等高線を描くと（関数を同じ高さのところで輪切りにすると）楕円になる.

この投稿では, 2変量正規分布の密度関数を同じ高さのところで輪切りにする方法を考える.

密度関数のとる値が同じになるのは,

$\displaystyle Z = (\boldsymbol{x}-\boldsymbol\mu)^{\top} \Lambda ({\boldsymbol x}-\boldsymbol \mu)$

の値が同じになることと同値である（他の部分は定数）.

分散共分散行列は実対称行列なので, 直交行列 $P$ を用いて,
$\displaystyle P^{\top} \Lambda P= \left( \begin{array}{c} \lambda_1 & 0\\ 0&\lambda_2 \\ \end{array} \right)$
と対角化できる（直交行列というのは転置行列が逆行列になる行列のこと）.

また,
$\displaystyle D= \left( \begin{array}{c} \lambda_1^{-1/2} & 0 \\ 0 & \lambda_2^{-1/2} \\ \end{array} \right)$
とすると,

$\displaystyle D P^{\top} \Lambda P D= \left( \begin{array}{c} 1 & 0 \\ 0 & 1 \\ \end{array} \right)$

である.

$\boldsymbol y = DP^{\top} (\boldsymbol {x} - \boldsymbol {\mu})$

と変換すると,

$\displaystyle Z = \boldsymbol{y}^{\top} \boldsymbol{y}$

である.

ベクトルの成分を,

$\displaystyle \boldsymbol y = \left( \begin{array}{c} y_1 \\ y_2 \\ \end{array} \right)$

として書き下すと,

$\displaystyle Z = y_1^2+y_2^2$

となり, これは円の方程式である.

このことは円の各点 $\boldsymbol y$ を $DP^{\top}$ で変換してやると楕円になることを意味する.

さて, いま考えているのは確率密度関数の等高線なので, どうせなら適当な確率 $p_0$ にもとづく高さで輪切りにしたい気持ちがある.

変換後の $\boldsymbol y$ の世界で考えると, $\boldsymbol y$ は独立な2つの標準正規分布（平均0, 分散1の正規分布）に従う確率変数である.

独立に標準正規分布に従うd個の確率変数をそれぞれ2乗して足し合わせた確率変数の従う分布を, 自由度dのカイ二乗分布という.

カイ二乗分布は統計学でよく使うのでたいていの統計ソフトでは $\Pr(Z \le z_0) = p_0$ を与える $z_0$ が計算できるようになっている.

$z_0$ の高さで密度関数を輪切りにしたものを確率楕円と呼ぶ.

実装例（Julia）

module Ellipse

using Distributions
using LinearAlgebra

function ellipse2d(μ, Λ, p, n)
    A = eigen(Λ)
    P = A.vectors
    D = diagm(sqrt.(inv.(A.values)))
    c = quantile(Chisq(2), p)
    Q = P'*D*sqrt(c)
    t = 0:(2pi/n):2pi
    x = zeros(length(t), 2)
    #z = zeros(length(t))
    for i in eachindex(t)
        x[i,:] = Q*[cos(t[i]), sin(t[i])] + μ #媒介変数表示
        #z[i] = (x[i,:]-μ)'*Λ*(x[i,:]-μ)
    end
    return x
end

#μなし（原点中心）
function ellipse2d(Λ, p, n)
    A = eigen(Λ)
    P = A.vectors
    D = diagm(sqrt.(inv.(A.values)))
    c = quantile(Chisq(2), p)
    Q = P'*D*sqrt(c)
    t = 0:(2pi/n):2pi
    x = zeros(length(t), 2)
    for i in eachindex(t)
        x[i,:] = Q*[cos(t[i]), sin(t[i])]
    end
    return x
end    

function countpoints(Y, x)
    p = zero(Int)
    for i in axes(Y,2)
        #どこかしら曲線より上でかつどこかしら曲線より下
        p += any(x[:,1] .< Y[1,i] .&& x[:,2] .< Y[2,i]) && any(Y[1,i] .< x[:,1] .&& Y[2,i] .< x[:,2])
    end
    return p    
end

end

適当に（乱数で）サンプルした2変量正規分布に従う点をプロットしてやると大体楕円におさまることがわかる.

ここでは $p_0$ に0.8を与えているが, 点の数を数えるとおおよそ期待値通りになっている.

julia> (m/20000)
0.79855

using Distributions
using Random
using LinearAlgebra
using Distributions
using Plots

rng = MersenneTwister(4126)
μ = randn(rng, 2)
Λ = Symmetric(rand(rng, Wishart(3, [1 0;0 1])))
X = rand(rng, MvNormal(μ, inv(Λ)), 20000)
x = Ellipse.ellipse2d(μ, Λ, 0.8, 201)

scatter(X[1,:], X[2,:], msw=0, ma=0.01, legend=false, tick_direction=:out)
plot!(x[:,1], x[:,2])
# png("./Desktop/scatter.png")

m = Ellipse.countpoints(X, x)
(m/20000)

補足：固有値分解

2変量正規分布の分散共分散行列は2つの分散 $\sigma_1^2, \sigma_2^2$ と相関係数 $\rho$ 用いて,

$\displaystyle \Sigma = \left( \begin{array}{c} \sigma_1^2 & \rho \sigma_1 \sigma_2\\ \rho \sigma_1 \sigma_2 &\lambda_2 \\ \end{array} \right)$

と表すこともできる.

$\Sigma$ の固有値を求めるため, 方程式
$\left| \begin{array}{c} 1 - \lambda & -\rho \\ -\rho & 1- \lambda \\ \end{array} \right| = 0$
を解くと,
$\lambda = 1 \pm \rho$
で, これが固有値である.

対応する固有ベクトルはそれぞれ,
$\displaystyle \frac{1}{\sqrt{2}}\left( \begin{array}{c} 1 \\ 1 \\ \end{array} \right)$
$\displaystyle \frac{1}{\sqrt{2}}\left( \begin{array}{c} -1 \\ 1 \\ \end{array} \right)$
であり, $\Sigma$ は
$P=\frac{1}{\sqrt{2}}\left( \begin{array}{c} 1 & -1\\ 1&1 \\ \end{array} \right)$
により
$P^\top \Sigma P= \left( \begin{array}{c} 1+\rho & 0 \\ 0&1-\rho \\ \end{array} \right)$
と対角化できる.