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) の推定を行うことをユニバーサルクリギングという。


  • 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を左辺に移動して (YXβ)=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}を次のように定義する。あるベクトルλ=(λ1,,λn)\lambda = \left(\lambda_{1} , \cdots ,\lambda_{n} \right)についてY(s0)Y \left( s_{0} \right)最良線形不偏予測量bLUP, Best Linear Unbiased Predictorλ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(x0XTΣ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

パート1. 条件付き期待値

E[(Y(s0)h(y))2y] 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))2y]=E[(Y(s0)E(Y(s0)y))2y]+[E(Y(s0)y)h(y)]2E[(Y(s0)E(Y(s0)y))2y] \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)である。


パート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)変量の正規分布に従うという仮定から始める。これをブロック行列形式で表すと次のようになる。 [Y0Y1Yn]=[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に対し多変量正規分布に従うランダムベクトルXNn(μ,Σ)\mathbf{X} \sim N_{n} \left( \mu , \Sigma \right)が与えられているとする。すると、条件付き確率ベクトルX1X2:ΩRm\mathbf{X}_{1} | \mathbf{X}_{2} : \Omega \to \mathbb{R}^{m}も依然として多変量正規分布に従い、具体的には次のように母平均ベクトルと母共分散行列を持つ。 X1X2Nm(μ1+Σ12Σ221(X2μ2),Σ11Σ12Σ221Σ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(YXβ)=x0Tβ+γTΣ1(YXβ) \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(x0XTΣ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. ↩︎