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라 한다. 이렇게 모델을 세워서 추정하는 행위 자체를 크리깅이라 부르기도 한다.

유니버설 크리깅

X:=[1X1(s1)⋯Xp(s1)1X1(s2)⋯Xp(s2)⋮⋮⋱⋮1X1(sn)⋯Xp(sn)] \mathbf{X} := \begin{bmatrix} 1 & X_{1} \left( s_{1} \right) & \cdots & X_{p} \left( s_{1} \right) \\ 1 & X_{1} \left( s_{2} \right) & \cdots & X_{p} \left( s_{2} \right) \\ \vdots & \vdots & \ddots & \vdots \\ 1 & X_{1} \left( s_{n} \right) & \cdots & X_{p} \left( s_{n} \right) \end{bmatrix} 독립변수 X1,⋯ ,XpX_{1} , \cdots , X_{p} 에 대해 위와 같은 디자인 매트릭스와 회귀계수의 벡터를 β=(β0,β1,⋯ ,βp)\beta = \left( \beta_{0} , \beta_{1} , \cdots , \beta_{p} \right) 라 할 때, 다변량정규분포를 따르는 ε∼Nn(0,Σ)\varepsilon \sim N_{n} \left( \mathbf{0} , \Sigma \right) 에 대해 다음의 모델 Y=Xβ+ε \mathbf{Y} = \mathbf{X} \beta + \varepsilon 을 통해 새로운 사이트site s0s_{0} 에서의 x(s0)=(x1(s0),⋯ ,xp(s0))\mathbf{x} \left( s_{0} \right) = \left( x_{1} \left( s_{0} \right) , \cdots , x_{p} \left( s_{0} \right) \right) 를 사용해 Y(s0)Y \left( s_{0} \right) 를 추정하는 과정을 유니버설 크리깅universal Kriging이라 한다.

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

설명

유니버설 크리깅이란 쉽게 말해 크리깅에 다중회귀분석이 가미된 크리깅을 말하며, 오디너리 크리깅에서의 μ\mu 는 유니버설 크리깅에서의 β0\beta_{0} 와 같은 역할을 한다. 시계열분석으로 치자면 동적 회귀 모델이라고 보면 되는데, 우리가 관심있는 변수인 YY 를 설명하기 위한 정보로써 순수하게 공간적인 구조 말고도 다른 변수들이 있을때 고려하게 된다.

수식을 보면 사실 회귀분석 파트가 좀 있어서 달라보일 뿐이지 Xβ\mathbf{X} \beta 을 좌변으로 넘겨서 (Y−Xβ)=0+ε \left( \mathbf{Y} - \mathbf{X} \beta \right) = \mathbf{0} + \varepsilon 이라 두면 그냥 평균벡터가 영벡터 μ=0\mu = \mathbf{0} 인 오디너리 크리깅 그 자체다.

공식

랜덤필드 {Y(sk)}k=1n\left\{ Y (s_{k}) \right\}_{k=1}^{n} 가 내재적 정상 공간과정이라고 하고 새롭계 예측하고 싶은 지점을 s0s_{0} 이라 하자. 배리오그램 2γ2 \gamma 에 대해 벡터 γ∈Rn\gamma \in \mathbb{R}^{n} 을 (γ)i:=(γ(s0−si))\left( \gamma \right)_{i} := \left( \gamma \left( s_{0} - s_{i} \right) \right) 와 같이 정의하자. 어떤 벡터 λ=(λ1,⋯ ,λn)\lambda = \left(\lambda_{1} , \cdots ,\lambda_{n} \right) 에 대해 Y(s0)Y \left( s_{0} \right) 의 최선선형불편예측량bLUP, Best Linear Unbiased Predictor은 λ\lambda 과 Y\mathbf{Y} 의 내적인 λTY=[λ1⋯λn][Y(s1)Y(s2)⋮Y(sn)]=∑k=1nλkY(sk) \lambda^{T} \mathbf{Y} = \begin{bmatrix} \lambda_{1} & \cdots & \lambda_{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} \lambda_{k} Y \left( s_{k} \right) 이며, 벡터 λ\lambda 는 구체적으로 다음과 같이 구해진다. λ=Σ−1γ+Σ−1X(XTΣ−1X)−1(x0−XTΣ−1γ) \lambda = \Sigma^{-1} \gamma + \Sigma^{-1} \mathbf{X} \left( \mathbf{X}^{T} \Sigma^{-1} \mathbf{X} \right)^{-1} \left( \mathbf{x}_{0} - \mathbf{X}^{T} \Sigma^{-1} \gamma \right)

유도 1

Part 1. 조건부 기대값

E[(Y(s0)−h(y))2∣y] E \left[ \left. \left( Y \left( s_{0} \right) - h \left( \mathbf{y} \right) \right)^{2} \right| \mathbf{y} \right] 우리의 목표는 데이터 y=(y(s1),⋯ ,y(sn))\mathbf{y} = \left( y \left( s_{1} \right) , \cdots , y \left( s_{n} \right) \right) 가 주어져 있을 때 위와 같은 에러를 최소화하는 h(y)h \left( \mathbf{y} \right) 를 구하는 것이다. E[(Y(s0)−h(y))2∣y]=E[(Y(s0)−E(Y(s0)∣y))2∣y]+[E(Y(s0)∣y)−h(y)]2≥E[(Y(s0)−E(Y(s0)∣y))2∣y] \begin{align*} & E \left[ \left. \left( Y \left( s_{0} \right) - h \left( \mathbf{y} \right) \right)^{2} \right| \mathbf{y} \right] \\ =& E \left[ \left. \left( Y \left( s_{0} \right) - E \left( Y \left( s_{0} \right) | \mathbf{y} \right) \right)^{2} \right| \mathbf{y} \right] + \left[ E \left( Y \left( s_{0} \right) | \mathbf{y} \right) - h \left( \mathbf{y} \right) \right]^{2} \\ \ge & E \left[ \left. \left( Y \left( s_{0} \right) - E \left( Y \left( s_{0} \right) | \mathbf{y} \right) \right)^{2} \right| \mathbf{y} \right] \end{align*} 이 식에서 등식이 성립하는 필요충분조건은 h(y)=E(Y(s0)∣y)h \left( \mathbf{y} \right) = E \left( Y \left( s_{0} \right) | \mathbf{y} \right) 이다.

Part 2. 조건부 다변량정규분포

기본적으로 주어진 데이터는 다변량정규분포 Nn(0,Σ)N_{n} \left( \mathbf{0} , \Sigma \right) 를 따른다는 전제 하에, 여기에 새로운 Y0:=Y(S0)∼N(μ0,σ02)Y_{0} := Y \left( S_{0} \right) \sim N \left( \mu_{0} , \sigma_{0}^{2} \right) 를 추가해도 여전히 (1+n)(1+n)변량 변량정규분포를 따른다는 가정에서 시작한다. 이를 블럭행렬꼴로 나타내보면 다음과 같다. [Y0Y1⋮Yn]=[Y0Y]∼N1+n([μ0μ],[σ02γTγΣ]) \begin{bmatrix} Y_{0} \\ Y_{1} \\ \vdots \\ Y_{n} \end{bmatrix} = \begin{bmatrix} Y_{0} \\ \mathbf{Y} \end{bmatrix} \sim N_{1+n} \left( \begin{bmatrix} \mu_{0} \\ \mu \end{bmatrix} , \begin{bmatrix} \sigma_{0}^{2} & \gamma^{T} \\ \gamma & \Sigma \end{bmatrix} \right)

다변량정규분포의 조건부 평균과 분산: X=[X1X2]:Ω→Rnμ=[μ1μ2]∈RnΣ=[Σ11Σ12Σ21Σ22]∈Rn×n \begin{align*} \mathbf{X} =& \begin{bmatrix} \mathbf{X}_{1} \\ \mathbf{X}_{2} \end{bmatrix} & : \Omega \to \mathbb{R}^{n} \\ \mu =& \begin{bmatrix} \mu_{1} \\ \mu_{2} \end{bmatrix} & \in \mathbb{R}^{n} \\ \Sigma =& \begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{bmatrix} & \in \mathbb{R}^{n \times n} \end{align*} 위와 같이 조던블럭폼으로 나타낸 X\mathbf{X}, μ\mu, Σ\Sigma 에 대해 다변량정규분포를 따르는 랜덤벡터 X∼Nn(μ,Σ)\mathbf{X} \sim N_{n} \left( \mu , \Sigma \right) 가 주어져 있다고 하자. 그러면 조건부확률벡터 X1∣X2:Ω→Rm\mathbf{X}_{1} | \mathbf{X}_{2} : \Omega \to \mathbb{R}^{m} 는 여전히 다변량정규분포를 따르며, 구체적으로 다음과 같이 모평균 벡터와 모공분산행렬을 가진다. X1∣X2∼Nm(μ1+Σ12Σ22−1(X2−μ2),Σ11−Σ12Σ22−1Σ21) \mathbf{X}_{1} | \mathbf{X}_{2} \sim N_{m} \left( \mu_{1} + \Sigma_{12} \Sigma_{22}^{-1} \left( \mathbf{X}_{2} - \mu_{2} \right) , \Sigma_{11} - \Sigma_{12} \Sigma_{22}^{-1} \Sigma_{21} \right)

따라서 데이터 Y\mathbf{Y} 가 주어졌을 때 Y0Y_{0} 의 기대값은 s0s_{0} 에서의 x0=x(s0)∈Rp\mathbf{x}_{0} = \mathbf{x} \left( s_{0} \right) \in \mathbb{R}^{p} 에 대해 E(Y(s0)∣Y)=μ0+γTΣ−1(Y−Xβ)=x0Tβ+γTΣ−1(Y−Xβ) \begin{align*} & E \left( Y \left( s_{0} \right) | \mathbf{Y} \right) \\ =& \mu_{0} + \gamma^{T} \Sigma^{-1} \left( \mathbf{Y} - \mathbf{X} \beta \right) \\ =& \mathbf{x}_{0}^{T} \beta + \gamma^{T} \Sigma^{-1} \left( \mathbf{Y} - \mathbf{X} \beta \right) \end{align*} 이다. 한편 회귀계수의 벡터 β\beta 는 다변량정규분포를 따르는 ε\varepsilon 의 각 성분이 독립이 아니라서 그냥 OLSordinary Least Square의 β^=(XTX)−1XTY\hat{\beta} = \left( \mathbf{X}^{T} \mathbf{X} \right)^{-1} \mathbf{X}^{T} \mathbf{Y} 를 쓸 수 없고 WLSweighted Least Square의 β^=(XTΣ−1X)−1XTΣ−1Y\hat{\beta} = \left( \mathbf{X}^{T} \Sigma^{-1} \mathbf{X} \right)^{-1} \mathbf{X}^{T} \Sigma^{-1} \mathbf{Y} 가 쓰여야한다. E(Y(s0)∣Y)=x0Tβ+γTΣ−1Y−γTΣ−1Xβ=x0Tβ+γTΣ−1Y−γTΣ−1Xβ=γTΣ−1Y+x0T(XTΣ−1X)−1XTΣ−1Y−γTΣ−1X(XTΣ−1X)−1XTΣ−1Y \begin{align*} & E \left( Y \left( s_{0} \right) | \mathbf{Y} \right) \\ =& \mathbf{x}_{0}^{T} \beta + \gamma^{T} \Sigma^{-1} \mathbf{Y} - \gamma^{T} \Sigma^{-1} \mathbf{X} \beta \\ =& \mathbf{x}_{0}^{T} \beta \\ & + \gamma^{T} \Sigma^{-1} \mathbf{Y} \\ & - \gamma^{T} \Sigma^{-1} \mathbf{X} \beta \\ =& \gamma^{T} \Sigma^{-1} \mathbf{Y} \\ & + \mathbf{x}_{0}^{T} \left( \mathbf{X}^{T} \Sigma^{-1} \mathbf{X} \right)^{-1} \mathbf{X}^{T} \Sigma^{-1} \mathbf{Y} \\ & - \gamma^{T} \Sigma^{-1} \mathbf{X} \left( \mathbf{X}^{T} \Sigma^{-1} \mathbf{X} \right)^{-1} \mathbf{X}^{T} \Sigma^{-1} \mathbf{Y} \end{align*} 여기서 Σ\Sigma 는 대칭행렬이므로 (Σ−1)T=Σ−1\left( \Sigma^{-1} \right)^{T} = \Sigma^{-1} 이고, Y\mathbf{Y} 로 묶어내면 다음과 같이 λTY\lambda^{T} \mathbf{Y} 꼴을 얻는다. γTΣ−1Y+x0T(XTΣ−1X)−1XTΣ−1Y−γTΣ−1X(XTΣ−1X)−1XTΣ−1Y=[Σ−1γ+Σ−1X(XTΣ−1X)−1(x0−XTΣ−1γ)]TY=λTY \begin{align*} & \gamma^{T} \Sigma^{-1} \mathbf{Y} \\ & + \mathbf{x}_{0}^{T} \left( \mathbf{X}^{T} \Sigma^{-1} \mathbf{X} \right)^{-1} \mathbf{X}^{T} \Sigma^{-1} \mathbf{Y} \\ & - \gamma^{T} \Sigma^{-1} \mathbf{X} \left( \mathbf{X}^{T} \Sigma^{-1} \mathbf{X} \right)^{-1} \mathbf{X}^{T} \Sigma^{-1} \mathbf{Y} \\ = & \left[ \Sigma^{-1} \gamma + \Sigma^{-1} \mathbf{X} \left( \mathbf{X}^{T} \Sigma^{-1} \mathbf{X} \right)^{-1} \left( \mathbf{x}_{0} - \mathbf{X}^{T} \Sigma^{-1} \gamma \right) \right]^{T} \mathbf{Y} \\ = & \lambda^{T} \mathbf{Y} \end{align*}

■

실제 분석에서 σ02\sigma_{0}^{2} 과 Σ\Sigma 의 추정치들은 세미배리오그램의 모형을 통해 다음과 같이 정해진다. σ02=σ2+τ2Σ=σ2H(ϕ)+τ2I \begin{align*} \sigma_{0}^{2} =& \sigma^{2} + \tau^{2} \\ \Sigma =& \sigma^{2} H \left( \phi \right) + \tau^{2} I \end{align*} 여기서 τ2\tau^{2} 는 너겟 이펙트 분산(이론과 달리 실제 데이터에서 거리에 상관없이 기본적으로 보게 되는 공분산성)고, II 는 항등행렬이다.

  1. Banerjee. (2015). Hierarchical Modeling and Analysis for Spatial Data(2nd Edition): p41~42. ↩︎