logo

공간데이터분석에서 크리깅이란? 📂통계적분석

공간데이터분석에서 크리깅이란?

모델

오디너리 크리깅

공간데이터분석에서 랜덤필드 Y=(Y(s1),⋯ ,Y(sn))\mathbf{Y} = \left( Y \left( s_{1} \right) , \cdots , Y \left( s_{n} \right) \right) 의 평균 μ∈R\mu \in \mathbb{R} 과 공분산행렬 Σ∈Rn×n\Sigma \in \mathbb{R}^{n \times n} 과 다변량정규분포를 따르는 ε∼Nn(0,Σ)\varepsilon \sim N_{n} \left( \mathbf{0} , \Sigma \right) 에 대해, 다음의 모델 Y=μ1+ε \mathbf{Y} = \mu \mathbf{1} + \varepsilon 을 통해 새로운 사이트site s0s_{0} 의 Y(s0)Y \left( s_{0} \right) 를 추정한 값을 오디너리 크리깅 추정치ordinary Kriging Estimate라 한다. 이렇게 모델을 세워서 추정하는 행위 자체를 크리깅이라 부르기도 한다.

  • 1=(1,⋯ ,1)\mathbf{1} = (1 , \cdots , 1) 은 모든 성분이 11 인 1벡터다.
  • Nn(0,Σ)N_{n} \left( \mathbf{0} , \Sigma \right) 는 다변량정규분포를 의미한다.

설명

어원

크리깅kringing이란 다니 G. 크리거daine G. Krige이라는 거장의 이름이 그대로 동사화 된 것이다. 보통의 통계학에서 예측, 예상, 적합, 추정 등등으로 쓰는 표현은 물론 ‘비어있는 공간’의 값을 채워넣는다는 의미에서 인터폴레이션 기법처럼 설명하는 경우도 제법 있는데, 이 모든 설명을 줄여서 크리깅하다라는 일반동사가 된 것으로 보면 된다.

좁은 의미에서의 응용수학, 컴퓨터 알고리즘, 머신러닝 기법 등과 구분되는 크리깅의 특징은 어쨌거나 통계학답게 그 평균(점추정량)만이 아니라 그 분산까지 고려한다는 것이다. 상상하기를, 각 지점에서 분산이 높은 곳의 크리깅 추정치는 마찬가지로 분산이 크고 분산이 낮은 곳들 사이의 지점에선 분산이 낮을 것이다. 이는 어떤 데이터의 관측소 입지를 선정하는데에도 쓰일 수 있는데, 예로써 미세먼지의 농도를 측정한다고 하면 미세먼지를 어떻게 측정할지가 궁금한 게 아니라 그 측정이 얼마나 정확한지―다시 말해 측정치의 분산이 가장 높은 곳을 선정하는 식의 접근법이 있다.

종속성

Y=μ1+ε,where ε∼Nn(0,Σ) \mathbf{Y} = \mu \mathbf{1} + \varepsilon \qquad , \text{where } \varepsilon \sim N_{n} \left( 0, \Sigma \right) 모델의 수식을 보면 회귀분석이나 시계열분석과 달리 ε\varepsilon 이야말로 우리의 관심사다. Σ\Sigma 가 대각행렬, 즉 관측치마다의 종속성이 없다면 애초에 공간적인 구조가 없다는 의미고 딱히 크리깅이라는 걸 할 이유가 없다. 실제 분석이서 이 Σ\Sigma 는 세미배리오그램의 모형을 통해 다음과 같이 정해진다. Σ=σ2H(ϕ)+τ2I \Sigma = \sigma^{2} H \left( \phi \right) + \tau^{2} I 여기서 τ2\tau^{2} 는 너겟 이펙트 분산(이론과 달리 실제 데이터에서 거리에 상관없이 기본적으로 보게 되는 공분산성)고, II 는 항등행렬이다.

일반화

다음과 같이 다른 독리변수에 대해 일반화 모델을 사용하는 크리깅을 유니버설 크리깅이라고 한다. Y=Xβ+ε \mathbf{Y} = \mathbf{X} \beta + \varepsilon

공식

랜덤필드 {Y(sk)}k=1n\left\{ Y (s_{k}) \right\}_{k=1}^{n} 가 내재적 정상 공간과정이라고 하고 새롭계 예측하고 싶은 지점을 s0s_{0} 이라 하자. 배리오그램 2γ2 \gamma 에 대해 행렬 Γ∈Rn×n\Gamma \in \mathbb{R}^{n \times n} 를 (Γ)ij:=γ(si−sj)\left( \Gamma \right)_{ij} := \gamma \left( s_{i} - s_{j} \right) 와 같이 정의하고 벡터 γ0∈Rn\gamma_{0} \in \mathbb{R}^{n} 을 (γ0)i:=(γ(s0−si))\left( \gamma_{0} \right)_{i} := \left( \gamma \left( s_{0} - s_{i} \right) \right) 와 같이 정의하자. 어떤 벡터 l=(l1,⋯ ,ln)l = \left( l_{1} , \cdots , l_{n} \right) 에 대해 Y(s0)Y \left( s_{0} \right) 의 최선선형불편예측량bLUP, Best Linear Unbiased Predictor은 ll 과 Y\mathbf{Y} 의 내적인 lTY=[l1⋯ln][Y(s1)Y(s2)⋮Y(sn)]=∑k=1nlkY(sk) l^{T} \mathbf{Y} = \begin{bmatrix} l_{1} & \cdots & l_{n} \end{bmatrix} \begin{bmatrix} Y \left( s_{1} \right)\\ Y \left( s_{2} \right)\\ \vdots\\ Y \left( s_{n} \right) \end{bmatrix} = \sum_{k=1}^{n} l_{k} Y \left( s_{k} \right) 이며, 벡터 ll 은 구체적으로 다음과 같이 구해진다. l=Γ−1(γ0+1−1TΓ−1γ01TΓ−111) l = \Gamma^{-1} \left( \gamma_{0} + {{ 1 - \mathbf{1}^T \Gamma^{-1} \gamma_{0} } \over { \mathbf{1}^T \Gamma^{-1} \mathbf{1} }} \mathbf{1} \right)

유도 1

구체적으로 크리깅이 어떻게 이루어지는지 수식으로 알아보는 과정이다. 이 공식의 가정에서 가우시안 프로세스라는 가정까지 추가되면 오디너리 크리깅이 된다.

Part 1. 최적화 문제

어떤 상수 l1,⋯ ,ln,δ0∈Rl_{1} , \cdots , l_{n} , \delta_{0} \in \mathbb{R} 들에 대해 새로운 Y(s0)Y \left( s_{0} \right) 를 기존 데이터의 선형결합 y^(s0)=l1y1+⋯+lnyn+δ0 \hat{y} \left( s_{0} \right) = l_{1} y_{1} + \cdots + l_{n} y_{n} + \delta_{0} 으로 예측하고 싶다고 하자. 이는 다시 말해 목적 함수 E[Y(s0)−(∑klkY(sk)+δ0)]2 E \left[ Y \left( s_{0} \right) - \left( \sum_{k} l_{k} Y \left( s_{k} \right) + \delta_{0} \right) \right]^{2} 를 최소화하는 최적해 l1,⋯ ,ln,δ0l_{1} , \cdots , l_{n} , \delta_{0} 를 찾는 것이다.

내재적 정상성의 정의: 유클리드 공간의 픽스된 부분집합 D⊂RrD \subset \mathbb{R}^{r} 에서 확률변수 Y(s):Ω→R1Y(s) : \Omega \to \mathbb{R}^{1} 의 집합인 공간과정 {Y(s)}s∈D\left\{ Y(s) \right\}_{s \in D} 와 방향벡터 h∈Rr\mathbf{h} \in \mathbb{R}^{r} 를 생각해보자. 구체적으로 n∈Nn \in \mathbb{N} 개의 사이트site를 {s1,⋯ ,sn}⊂D\left\{ s_{1} , \cdots , s_{n} \right\} \subset D 과 같이 나타내고, Y(s)Y(s) 는 모든 s∈Ds \in D 에 대해 분산이 존재하는 것으로 가정한다. [Y(s+h)−Y(s)]\left[ Y \left( s + \mathbf{h} \right) - Y(s) \right] 의 평균이 00 이면서 분산이 오직 h\mathbf{h} 에만 의존하면 {Y(sk)}\left\{ Y \left( s_{k} \right) \right\} 가 내재적 정상성intrinsic Stationarity을 가진다고 한다. E[Y(s+h)−Y(s)]=0Var⁡[Y(s+h)−Y(s)]=2γ(h) \begin{align*} E \left[ Y \left( s + \mathbf{h} \right) - Y(s) \right] =& 0 \\ \Var \left[ Y \left( s + \mathbf{h} \right) - Y(s) \right] =& 2 \gamma ( \mathbf{h} ) \end{align*}

여기서 {Y(sk)}k=1n\left\{ Y \left( s_{k} \right) \right\}_{k=1}^{n} 가 내재적 정상성을 가진다면, ∑klk=1\sum_{k} l_{k} = 1 이라 제약조건을 둠으로써 E[Y(s0)−(∑klkY(sk))]=E[∑klkY(s0)−(∑klkY(sk))]=∑klkE[Y(s0)−Y(sk)]=0 \begin{align*} & E \left[ Y \left( s_{0} \right) - \left( \sum_{k} l_{k} Y \left( s_{k} \right) \right) \right] \\ =& E \left[ \sum_{k} l_{k} Y \left( s_{0} \right) - \left( \sum_{k} l_{k} Y \left( s_{k} \right) \right) \right] \\ =& \sum_{k} l_{k} E \left[ Y \left( s_{0} \right) - Y \left( s_{k} \right) \right] \\ =& 0 \end{align*} 이 되게끔 할 수 있다. 이에 따라 우리가 최소화할 목적 함수는 δ0\delta_{0} 만 밖으로 빠진 다음의 꼴이 된다. E[Y(s0)−∑klkY(sk)]2+δ02 E \left[ Y \left( s_{0} \right) - \sum_{k} l_{k} Y \left( s_{k} \right) \right]^{2} + \delta_{0}^{2} 여기서 δ02\delta_{0}^{2} 는 예측과 상관없다. 실제로도 모델이 Y=μ1+ε\mathbf{Y} = \mu \mathbf{1} + \varepsilon 이라면 δ0\delta_{0} 는 μ\mu 에 해당해서 δ0=0\delta_{0} = 0 으로 두어 E[Y(s0)−∑klkY(sk)]2 E \left[ Y \left( s_{0} \right) - \sum_{k} l_{k} Y \left( s_{k} \right) \right]^{2} 라 해도 무방하다. 이제 a0=1a_{0} = 1 이라 두고 ak=−lka_{k} = - l_{k} 라 하면 E[Y(s0)−∑k=1nlkY(sk)]2=E[∑k=0nakY(sk)]2 E \left[ Y \left( s_{0} \right) - \sum_{k=1}^{n} l_{k} Y \left( s_{k} \right) \right]^{2} = E \left[ \sum_{k=0}^{n} a_{k} Y \left( s_{k} \right) \right]^{2} 이므로, 우리는 다음과 같은 최적화문제를 라그랑주 승수법으로 풀게 된 것이다. MinimizeE[∑k=0nakY(sk)]2subject to∑k=0nak=0 \begin{matrix} \text{Minimize} & \displaystyle E \left[ \sum_{k=0}^{n} a_{k} Y \left( s_{k} \right) \right]^{2} \\ \text{subject to} & \displaystyle \sum_{k=0}^{n} a_{k} = 0 \end{matrix}

Part 2. 세미배리오그램 γ\gamma

이제 E[∑k=0nakY(sk)]2E \left[ \sum_{k=0}^{n} a_{k} Y \left( s_{k} \right) \right]^{2} 를 ‘우리가 데이터에서 계산할 수 있을 것으로 가정할 수 있는’ 세미배리오그램에 종속된 꼴로 나타내보자.

세미배리오그램의 정의: 유클리드 공간의 픽스된 부분집합 D⊂RrD \subset \mathbb{R}^{r} 에서 확률변수 Y(s):Ω→R1Y(s) : \Omega \to \mathbb{R}^{1} 의 집합인 공간과정 {Y(s)}s∈D\left\{ Y(s) \right\}_{s \in D} 와 방향벡터 h∈Rr\mathbf{h} \in \mathbb{R}^{r} 를 생각해보자. 구체적으로 n∈Nn \in \mathbb{N} 개의 사이트를 {s1,⋯ ,sn}⊂D\left\{ s_{1} , \cdots , s_{n} \right\} \subset D 과 같이 나타내고, Y(s)Y(s) 는 모든 s∈Ds \in D 에 대해 분산이 존재하는 것으로 가정한다. 다음과 같이 정의되는 2γ(h)2 \gamma ( \mathbf{h} ) 를 배리오그램variogram이라 한다. 2γ(h):=E[Y(s+h)−Y(s)]2 2 \gamma ( \mathbf{h} ) := E \left[ Y \left( s + \mathbf{h} \right) - Y(s) \right]^{2} 특히 배리오그램의 절반 γ(h)\gamma ( \mathbf{h} ) 를 세미배리오그램semivariogram이라 한다.

γ(si−sj)\gamma \left( s_{i} - s_{j} \right) 를 두 사이트 si,sjs_{i}, s_{j} 사이의 방향벡터에 따른 세미배리오그램이라 하자. ∑0=1nak=0\sum_{0=1}^{n} a_{k} = 0 을 만족하는 아무 집합 {ak:k=1,⋯ ,n}⊂R\left\{ a_{k} : k = 1 , \cdots , n \right\} \subset \mathbb{R} 에 대해 다음이 성립한다. ∑i∑jaiajγ(si−sj)=−E[∑iaiY(si)]2 \sum_{i} \sum_{j} a_{i} a_{j} \gamma \left( s_{i} - s_{j} \right) = - E \left[ \sum_{i} a_{i} Y \left( s_{i} \right) \right]^{2} 이는 다음의 전개 ∑i∑jaiajγ(si−sj)=12∑i∑jaiajVar⁡[Y(si)−Y(sj)]=12∑i∑jaiajE[Y(si)−Y(sj)]2=12∑i∑jaiajE([Y(si)]2−2Y(si)Y(sj)+[Y(sj)]2)=−∑i∑jaiajE(Y(si)Y(sj))∵cases of i=j=−E∑i∑jaiajY(si)Y(sj)=−E∑iaiY(si)∑jajY(sj)=−E[∑iaiY(si)]2 \begin{align*} & \sum_{i} \sum_{j} a_{i} a_{j} \gamma \left( s_{i} - s_{j} \right) \\ =& {{ 1 } \over { 2 }} \sum_{i} \sum_{j} a_{i} a_{j} \Var \left[ Y \left( s_{i} \right) - Y \left( s_{j} \right) \right] \\ =& {{ 1 } \over { 2 }} \sum_{i} \sum_{j} a_{i} a_{j} E \left[ Y \left( s_{i} \right) - Y \left( s_{j} \right) \right]^{2} \\ =& {{ 1 } \over { 2 }} \sum_{i} \sum_{j} a_{i} a_{j} E \left( \left[ Y \left( s_{i} \right) \right]^{2} - 2 Y \left( s_{i} \right) Y \left( s_{j} \right) + \left[ Y \left( s_{j} \right) \right]^{2} \right) \\ =& - \sum_{i} \sum_{j} a_{i} a_{j} E \left( Y \left( s_{i} \right) Y \left( s_{j} \right) \right) & \because \text{cases of } i = j \\ =& - E \sum_{i} \sum_{j} a_{i} a_{j} Y \left( s_{i} \right) Y \left( s_{j} \right) \\ =& - E \sum_{i} a_{i} Y \left( s_{i} \right) \sum_{j} a_{j} Y \left( s_{j} \right) \\ =& - E \left[ \sum_{i} a_{i} Y \left( s_{i} \right) \right]^{2} \end{align*} 에서 확인할 수 있다. 이제 γij=γ(si−sj)\gamma_{ij} = \gamma \left( s_{i} - s_{j} \right) 라 하고 γ0j=γ(s0−sj)\gamma_{0j} = \gamma \left( s_{0} - s_{j} \right) 라 두면 우리의 목적함수는 E[∑iaiY(si)]2=−∑i∑jaiajγ(si−sj)=−∑i∑jliljγij+2∑iliγ0i \begin{align*} & E \left[ \sum_{i} a_{i} Y \left( s_{i} \right) \right]^{2} \\ =& - \sum_{i} \sum_{j} a_{i} a_{j} \gamma \left( s_{i} - s_{j} \right) \\ =& - \sum_{i} \sum_{j} l_{i} l_{j} \gamma_{ij} + 2 \sum_{i} l_{i} \gamma_{0i} \end{align*} 와 같이 나타나고, 라그랑주 승수법에 따라 제약조건 ∑ili=1\sum_{i} l_{i} = 1 에 라그랑주 승수lagrange Multiplier를 곱해서 빼면  −∑i∑jliljγij+2∑iliγ0i−λ∑ili \ - \sum_{i} \sum_{j} l_{i} l_{j} \gamma_{ij} + 2 \sum_{i} l_{i} \gamma_{0i} - \lambda \sum_{i} l_{i} 을 얻는다. 따라서 lil_{i} 들로 편미분해서  −∑jljγij+γ0i−λ=0 \ - \sum_{j} l_{j} \gamma_{ij} + \gamma_{0i} - \lambda = 0 이 될 때 목적함수를 최소화 시킬 수 있음을 알 수 있다.

Part3. 최적해

이제 구체적인 최적해의 공식을 유도할 것이다. 행렬 Γ∈Rn×n\Gamma \in \mathbb{R}^{n \times n} 의 (i,j)(i,j) 성분을 γij\gamma_{ij} 로, 즉 (Γ)ij:=γij\left( \Gamma \right)_{ij} := \gamma_{ij} 라 하고 벡터 γ0∈Rn\gamma_{0} \in \mathbb{R}^{n} 을 (γ0)i:=γ0i\left( \gamma_{0} \right)_{i} := \gamma_{0i} 와 같이 정의하자. 계수의 벡터 역시 l:=(l1,⋯ ,ln)∈Rnl := \left( l_{1} , \cdots , l_{n} \right) \in \mathbb{R}^{n} 이라 두면 Part2에서 얻은 식은 다음과 같은 행렬/벡터 폼으로 나타낼 수 있다. −∑jljγij+γ0i−λ=0  ⟹  −Γl+γ0−λ1=0  ⟹  Γl+λ1=γ0 \begin{align*} & - \sum_{j} l_{j} \gamma_{ij} + \gamma_{0i} - \lambda = 0 \\ \implies & - \Gamma l + \gamma_{0} - \lambda \mathbf{1} = 0 \\ \implies & \Gamma l + \lambda \mathbf{1} = \gamma_{0} \end{align*} 한편 제약조건에서 ∑ili=1  ⟺  [1⋯1][l1⋮ln]=1  ⟺  1Tl=1 \sum_{i} l_{i} = 1 \iff \begin{bmatrix} 1 & \cdots & 1 \end{bmatrix} \begin{bmatrix} l_{1} \\ \vdots \\ l_{n} \end{bmatrix} = 1 \iff \mathbf{1}^{T} l = 1 이므로 1Tl=1\mathbf{1}^{T} l = 1 도 얻을 수 있다. 여기서 xT\mathbf{x}^{T} 는 x\mathbf{x} 의 트랜스포즈를 나타낸다. 우선 단독으로 λ\lambda 는 1=1Tl=1TΓ−1Γl=1TΓ−1(γ0−λ1)=1TΓ−1γ0−λ1TΓ−11 \begin{align*} 1 =& \mathbf{1}^T l \\ =& \mathbf{1}^T \Gamma^{-1} \Gamma l \\ =& \mathbf{1}^T \Gamma^{-1} \left( \gamma_{0} - \lambda \mathbf{1} \right) \\ =& \mathbf{1}^T \Gamma^{-1} \gamma_{0} - \lambda \mathbf{1}^T \Gamma^{-1} \mathbf{1} \end{align*} 이고, 정리하면  −λ=1−1TΓ−1γ01TΓ−11 \ - \lambda = {{ 1 - \mathbf{1}^T \Gamma^{-1} \gamma_{0} } \over { \mathbf{1}^T \Gamma^{-1} \mathbf{1} }} 와 같이 나타낼 수 있다. 이제 ll 은 다 구한 것이나 마찬가지다. Γl+λ1=γ0  ⟹  Γl=γ0−λ1  ⟹  Γl=γ0+1−1TΓ−1γ01TΓ−111  ⟹  l=Γ−1(γ0+1−1TΓ−1γ01TΓ−111) \begin{align*} & \Gamma l + \lambda \mathbf{1} = \gamma_{0} \\ \implies & \Gamma l = \gamma_{0} - \lambda \mathbf{1} \\ \implies & \Gamma l = \gamma_{0} + {{ 1 - \mathbf{1}^T \Gamma^{-1} \gamma_{0} } \over { \mathbf{1}^T \Gamma^{-1} \mathbf{1} }} \mathbf{1} \\ \implies & l = \Gamma^{-1} \left( \gamma_{0} + {{ 1 - \mathbf{1}^T \Gamma^{-1} \gamma_{0} } \over { \mathbf{1}^T \Gamma^{-1} \mathbf{1} }} \mathbf{1} \right) \end{align*}

■

  1. Banerjee. (2015). Hierarchical Modeling and Analysis for Spatial Data(2nd Edition): p25, 40~41. ↩︎