11 Sep 2017 » probability, statistics
2017-09-11-probability_1

수식이 깨질 경우 새로고침을 눌러주세요.

Probability for Data-Science-13

다변수 확률 모형 ②

  • Categorical (Multinoulli)
  • Multinomial
  • ✓Dirichlet
  • Multivariate Gaussian Normal

Dirichlet distribution (디리클레 분포)

디리클레 분포(Dirichlet distribution)는 베타 분포의 다차원 확장

  • 베타 분포 → 0과 1사이의 값을 가지는 단일(univariate) 확률 변수의 베이지안 모형에 사용된다.
  • 디리클레 분포는 0과 1사이의 사이의 값을 가지는 다변수(multivariate) 확률 변수의 베이지안 모형에 사용된다.

다만, 디리클레 분포는 다변수 확률 변수들의 합이 1이되어야 한다는 제한 조건을 가진다.

참고 : 베타분포

디리클레 분포의 확률 밀도 함수 $$ f(x_1, x_2, \cdots, x_K) = \frac{1}{\mathrm{B}(\boldsymbol\alpha)} \prod_{i=1}^K x_i^{\alpha_i - 1} $$

$$ \begin{eqnarray} 단,\;\; \mathrm{B}(\boldsymbol\alpha) &=& \frac{\prod_{i=1}^K \Gamma(\alpha_i)} {\Gamma\bigl(\sum_{i=1}^K \alpha_i\bigr)} \\ \\&=& \frac{\Gamma(\alpha_1)\cdots\Gamma(\alpha_K)}{\Gamma(\alpha_1+\cdots+\alpha_K)} \end{eqnarray} $$$$또한,\;\;\; (0 \leq x_i \leq 1)\;이며\;\; \sum_{i=1}^{K} x_i = 1 $$$$ 이\;\; 식에서 \;\; \boldsymbol\alpha =\; (\alpha_1, \alpha_2, \ldots, \alpha_K)는 \;\;디리클레 \;\;분포의 \;\;모수 \;\;벡터이다. $$

ex) K=3인 디리클레 분포를 따르는 확률 변수는 다음과 같은 값을 샘플로 가질 수 있다.

$$(1, 0, 0), \;\;\;(0.5, 0.5, 0),\;\;\; (0.2, 0.3, 0.5)$$

베타 분포와 디리클레 분포의 관계

베타 분포는 $K=2$ 인 디리클레 분포라고 볼 수 있다.

$$ f(x_1, x_2) = \frac{1}{\mathrm{B}(\boldsymbol\alpha)} \prod_{i=1}^2 x_i^{\alpha_i - 1} $$
$$ \begin{eqnarray} f(x_1, x_2) &=& \frac{1}{\mathrm{B}(\boldsymbol\alpha)} \prod_{i=1}^2 x_i^{\alpha_i - 1} \\ &=& \frac{1}{\mathrm{B}(\alpha_1,\alpha_2)} x_1^{\alpha_1 - 1}x_2^{\alpha_2 - 1} \end{eqnarray} $$
$$ 여기서\;\; x_2 \;=\; 1 - x_1\; 이고, \alpha_1= a,\;\; \alpha_2= b\; 라 \;\;하면 $$
$$ \begin{eqnarray} \frac{1}{\mathrm{B}(\alpha_1,\alpha_2)} x_1^{\alpha_1 - 1}x_2^{\alpha_2 - 1} &=& \frac{1}{\mathrm{B}(a,b)} x_1^{a - 1}(1-x_1)^{b - 1} \\ \\&=& \text{Beta}(x_1;a,b) \end{eqnarray} $$

디리클레 분포의 모멘트 특성

기대값 $$ \begin{eqnarray} E[x_k] &=& \dfrac{\alpha_k}{\boldsymbol\alpha} \\ &=& \dfrac{\alpha_k}{\sum\alpha_k} \end{eqnarray} $$

$$ 위에서\;\; 보듯이\;\; \boldsymbol{\alpha} = \sum\alpha_k $$$$ 기대값\;\; 공식을\;\; 보면\;\; 모수인\;\; \boldsymbol\alpha \;=\; (\alpha_1,\; \alpha_2,\; \ldots,\; \alpha_K)는\\ (x_1, x_2, \ldots, x_K)\;\; 중\;\; 어느\;\; 수가\;\; 더\;\; 크게\;\; 나올\;\; 가능성이\;\; 높은\;지를\;\; 결정하는\;\; 형상\;\; 인자(shape \;factor)임을\;\; 알\;\; 수\;\; 있다.\;\;\\ 모든\;\; \alpha_i값이\;\; 동일하면\;\; 모든\;\; x_i의\;\; 분포가\;\; 같아진다. $$

모드

$$ \dfrac{\alpha_k - 1}{\boldsymbol\alpha - K}$$

분산

$$\text{Var}[x_k] =\dfrac{\alpha_k(\boldsymbol\alpha - \alpha_k)}{\boldsymbol\alpha^2(\boldsymbol\alpha + 1)}$$
$$ 또한\;\; 분산\;\; 공식을\;\; 보면\;\; \boldsymbol\alpha의\;\; 절대값이\;\; 클수록\;\; 분산이\;\; 작아진다.\\ 즉,\;\; 어떤\;\; 특정한\;\; 값이\;\; 나올\;\; 가능성이\;\; 높아진다. $$

디리클레 분포의 응용

$$ 다음과\;\; 같은\;\; 문제를\;\; 보자.\;\; 이\;\; 문제는\;\; K=3이고\;\; \alpha_1 = \alpha_2 = \alpha_3\; 인\;\; Dirichlet\;\; 분포의\;\; 특수한\;\; 경우이다. $$
$$ x,\; y,\; z가\;\; 난수일\;\; 때,\;\; 항상\; x\;+\;y\;+\;z \;=\; 1이\;\; 되도록\;\; 하려면?\;\; (단,\;\; x,\;y,\;z는\;\; 균등하게\;\; 뽑아야\;\; 하며,\;\; x>0\;, y>0\;, z>0) $$

Sol

3차원 디리클레 문제는 3차원 공간 상에서 (1,0,0), (0,1,0), (0,0,1) 세 점을 연결하는 정삼각형 면위의 점을 생성하는 문제라고 볼 수 있다.

In [1]:
# -*- coding: utf-8 -*-

import seaborn as sns
import pandas as pd
import scipy as sp
import matplotlib as mpl
import matplotlib.pylab as plt
import numpy as np
%matplotlib inline
In [2]:
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

# 3차원 도화지를 꺼내라.
fig = plt.figure(figsize=(8,5))
ax = Axes3D(fig)


x = [1,0,0]
y = [0,1,0]
z = [0,0,1]

# x,y,z벡터의 각 원소를 zip으로 묶어서, 3차 좌표로 변환
## zip으로 좌표를 묶은 후 리스트가 두개가 있어야 그래프가 그려진다. 차원을 맞춰 주어야한다.
verts = [list(zip(x, y, z))]
ax.add_collection3d(Poly3DCollection(verts, edgecolor="k", lw=5, alpha=0.4))

ax.text(1, 0, 0, "(1,0,0)", position=(0.7,0.1))
ax.text(0, 1, 0, "(0,1,0)", position=(0,1.04))
ax.text(0, 0, 1, "(0,0,1)", position=(-0.2,0))

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

ax.set_xticks([])
ax.set_yticks([])
ax.set_zticks([])

ax.view_init(30, -20)
tmp_planes = ax.zaxis._PLANES 
# set origin ( http://stackoverflow.com/questions/15042129/changing-position-of-vertical-z-axis-of-3d-plot-matplotlib )
ax.yaxis._PLANES = (tmp_planes[2], tmp_planes[3], 
                    tmp_planes[0], tmp_planes[1], 
                    tmp_planes[4], tmp_planes[5])
ax.zaxis._PLANES = (tmp_planes[2], tmp_planes[3], 
                    tmp_planes[0], tmp_planes[1], 
                    tmp_planes[4], tmp_planes[5])
plt.show()

주의!

In [3]:
x = [1,0,0]
y = [0,1,0]
z = [0,0,1]

verts_1 = [zip(x, y, z)]
verts_2 = list(zip(x, y, z))
In [4]:
verts_1
Out[4]:
[<zip at 0x10f071288>]
In [5]:
verts_2
Out[5]:
[(1, 0, 0), (0, 1, 0), (0, 0, 1)]

위에서 보듯이 리스트와 zip은 같이 다녀야 한다!

In [6]:
x = [1,0,0]
y = [0,1,0]
z = [0,0,1]

# x,y,z벡터의 각 원소를 zip으로 묶어서, 3차 좌표로 변환
## zip으로 좌표를 묶은 후 리스트가 두개가 있어야 그래프가 그려진다. 차원을 맞춰 주어야한다.
verts_3 = [list(zip(x, y, z))]
In [7]:
verts_3
Out[7]:
[[(1, 0, 0), (0, 1, 0), (0, 0, 1)]]

그래프를 그릴 때, 위 처럼 리스트가 2개가 있어야한다. 차원을 맞춰줘야 한다!

다음 함수는 생성된 점들을 2차원 삼각형 위에서 볼 수 있도록 그려주는 함수이다.

In [8]:
def plot_triangle(X, kind):
    n1 = np.array([1, 0, 0])
    n2 = np.array([0, 1, 0])
    n3 = np.array([0, 0, 1])
    
    # n1과 n2의 중점
    n12 = (n1 + n2)/2
    m1 = np.array([1, -1, 0])
    m2 = n3 - n12
    
    # 노말라이징해서 단위벡터로 변환
    m1 = m1/np.linalg.norm(m1)
    m2 = m2/np.linalg.norm(m2)

    X1 = (X-n12).dot(m1)
    X2 = (X-n12).dot(m2)
    
    g = sns.jointplot(X1, X2, kind=kind, xlim=(-0.8,0.8), ylim=(-0.45,0.9))
    g.ax_joint.axis("equal")
    plt.show()

만약 이 문제를 단순하게 생각하여 서로 독립인 0과 1사이의 유니폼 확률 변수를 3개 생성하고 이들의 합이 1이 되도록 크기를 정규화(normalize)하면

다음 그림과 같이 삼각형의 중앙 근처에 많은 확률 분포가 집중된다. 즉, 확률 변수가 골고루 분포되지 않는다.

In [9]:
X1 = np.random.rand(1000, 3)
X1 = X1/X1.sum(axis=1)[:, np.newaxis]
plot_triangle(X1, kind="scatter")
In [10]:
plot_triangle(X1, kind="hex")

그러나 alpha=(1,1,1)인 디리클레 분포는 다음과 같이 골고루 샘플을 생성한다.

In [11]:
X2 = sp.stats.dirichlet((1,1,1)).rvs(1000)
plot_triangle(X2, kind="scatter")
In [12]:
plot_triangle(X2, kind="hex")

alpha가 (1,1,1) 이 아닌 경우

다음과 같이 특정 위치에 분포가 집중되도록 할 수 있다.

이 특성을 이용하면 다항 분포의 모수를 추정하는 베이지안 추정 문제에 응용할 수 있다.

In [13]:
def project(x):
    n1 = np.array([1, 0, 0])
    n2 = np.array([0, 1, 0])
    n3 = np.array([0, 0, 1])
    n12 = (n1 + n2)/2
    m1 = np.array([1, -1, 0])
    m2 = n3 - n12
    m1 = m1/np.linalg.norm(m1)
    m2 = m2/np.linalg.norm(m2)
    return np.dstack([(x-n12).dot(m1), (x-n12).dot(m2)])[0]

def project_reverse(x):
    n1 = np.array([1, 0, 0])
    n2 = np.array([0, 1, 0])
    n3 = np.array([0, 0, 1])
    n12 = (n1 + n2)/2
    m1 = np.array([1, -1, 0])
    m2 = n3 - n12
    m1 = m1/np.linalg.norm(m1)
    m2 = m2/np.linalg.norm(m2)
    return x[:,0][:, np.newaxis] * m1 + x[:,1][:, np.newaxis] * m2 + n12

eps = np.finfo(float).eps * 10
X = project([[1-eps,0,0], [0,1-eps,0], [0,0,1-eps]])
In [14]:
import matplotlib.tri as mtri
triang = mtri.Triangulation(X[:,0], X[:,1], [[0, 1, 2]])
refiner = mtri.UniformTriRefiner(triang)
triang2 = refiner.refine_triangulation(subdiv=6)
XYZ = project_reverse(np.dstack([triang2.x, triang2.y, 1-triang2.x-triang2.y])[0])
In [15]:
# alpha = (1,1,1) 일 때
pdf = sp.stats.dirichlet((1,1,1)).pdf(XYZ.T)

# 삼각 등고선을 그려라.
plt.tricontourf(triang2, pdf)
plt.axis("equal")
plt.show()
In [16]:
# (3, 4, 2) 일 때
pdf = sp.stats.dirichlet((3,4,2)).pdf(XYZ.T)
plt.tricontourf(triang2, pdf)
plt.axis("equal")
plt.show()
In [17]:
# (16, 24, 14) 일 때
pdf = sp.stats.dirichlet((16,24,14)).pdf(XYZ.T)
plt.tricontourf(triang2, pdf)
plt.axis("equal")
plt.show()

alpha 값에 따라 자료의 집중도가 달라진다.

모수의 추정

$$ \hat{\theta}\; : \;\;\; Hyper\;\; Parameter\;라고도\;\; 하며,\;\; 분포의\;\; 파라미터를\;\; 추정하는\;\; 파라미터이다. $$

베타분포

  • a = b = 1
  • a, b를 조작해서 어느 부분의 값을 많이 나오게 할 수 있다.

디리클레 분포 $$ x_1 + x_2 + x_3 = K \rightarrow x_1,\;x_2가\;\; 정해지면\;\; 나머지\;\; x_3는\;\; 자동으로 \;\;정해진다. $$

Parameter가 많을 수록 비모수적 방법을 사용한다.

  • seaborn histogram을 사용.

    자료의 분포를 히스토그램으로 구간별로 쌓아서 그리도록 한다. 앞에서 해왔던 시뮬레이션의 방법과 동일하다.

모수적 방법(Parameteric) VS 비모수적 방법(Non-Parameteric)

  • 일반적으로 파라미터가 5개, 10개가 넘어가면 비모수라고 한다.

Related Posts