📂統計的分析

空間データ分析におけるクリギングとは？

モデル

オーディナリークリギング

空間データ分析で、ランダムフィールド $\mathbf{Y} = \left( Y \left( s_{1} \right) , \cdots , Y \left( s_{n} \right) \right)$ の平均 $\mu \in \mathbb{R}$ と共分散行列 $\Sigma \in \mathbb{R}^{n \times n}$ が多変量正規分布に従う $\varepsilon \sim N_{n} \left( \mathbf{0} , \Sigma \right)$ ことについて、次のモデル $$ \mathbf{Y} = \mu \mathbf{1} + \varepsilon $$ を通じて、新しいサイト^site $s_{0}$ の $Y \left( s_{0} \right)$ を推定した値をオーディナリークリギング推定値^{ordinary Kriging Estimate}と呼ぶ。このようにモデルを立てて推定する行為自体をクリギングとも呼ぶ。

---

• $\mathbf{1} = (1 , \cdots , 1)$ は全ての成分が $1$ の1ベクトルだ。
• $N_{n} \left( \mathbf{0} , \Sigma \right)$ は多変量正規分布を意味する。

説明

語源

クリギング^krigingは、ダニエル・G・クリージ^{daine G. Krige}という巨匠の名前がそのまま動詞化されたものである。通常の統計学で予測、予想、適合、推定などと使う表現はもちろん、「空白の空間」の値を埋めるという意味で補間技術のように説明する場合もかなりあるが、これら全ての説明を短くしてクリギングするという一般動詞になったと考えられる。

狭義の応用数学、コンピュータアルゴリズム、機械学習技術などと区別されるクリギングの特徴は、とにかく統計学らしくその平均（点推定量）だけでなくその分散まで考慮することである。想像してみてほしいが、各地点で分散が高い場所のクリギング推定値は同様に分散が大きく、分散が低い場所同士の地点では分散が低いだろう。これはあるデータの観測所の位置を選定するのにも使われるが、例えば微粒子の濃度を測定するとしたら、微粒子をどのように測定するかが気になるのではなく、その測定がどれだけ正確か―つまり、測定値の分散が最も高い場所を選定するようなアプローチがある。

依存性

$$ \mathbf{Y} = \mu \mathbf{1} + \varepsilon \qquad , \text{where } \varepsilon \sim N_{n} \left( 0, \Sigma \right) $$ モデルの数式を見ると、回帰分析や時系列分析と異なり $\varepsilon$ こそが我々の関心事である。$\Sigma$ が対角行列、つまり観測値ごとの依存性がなければ、そもそも空間的な構造がないという意味であり、わざわざクリギングをする理由がない。実際の分析ではこの $\Sigma$ はセミバリオグラムのモデルを通じて次のように決定される。 $$ \Sigma = \sigma^{2} H \left( \phi \right) + \tau^{2} I $$ ここで $\tau^{2}$ はナゲット効果分散（理論と異なり実際のデータで距離に関係なく基本的に見られる共分散性）であり、$I$ は単位行列である。

一般化

次のように他の独立変数に対して一般化モデルを使用するクリギングをユニバーサルクリギングと呼ぶ。 $$ \mathbf{Y} = \mathbf{X} \beta + \varepsilon $$

公式

ランダムフィールド $\left\{ Y (s_{k}) \right\}_{k=1}^{n}$ が固有的定常空間過程であるとし、新たに予測したい地点を $s_{0}$ とする。バリオグラム $2 \gamma$ に対して行列 $\Gamma \in \mathbb{R}^{n \times n}$ を $\left( \Gamma \right)_{ij} := \gamma \left( s_{i} - s_{j} \right)$ のように定義し、ベクトル $\gamma_{0} \in \mathbb{R}^{n}$ を $\left( \gamma_{0} \right)_{i} := \left( \gamma \left( s_{0} - s_{i} \right) \right)$ のように定義する。あるベクトル $l = \left( l_{1} , \cdots , l_{n} \right)$ に対して $Y \left( s_{0} \right)$ の最良線形不偏予測量^{bLUP, Best Linear Unbiased Predictor}は $l$ と $\mathbf{Y}$ の内積であり、 $$ 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) $$ ベクトル $l$ は具体的に次のように求められる。 $$ 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. 最適化問題

ある定数 $l_{1} , \cdots , l_{n} , \delta_{0} \in \mathbb{R}$ に対して新しい $Y \left( s_{0} \right)$ を既存データの線形結合 $$ \hat{y} \left( s_{0} \right) = l_{1} y_{1} + \cdots + l_{n} y_{n} + \delta_{0} $$ で予測したいとする。これはつまり、目的関数 $$ E \left[ Y \left( s_{0} \right) - \left( \sum_{k} l_{k} Y \left( s_{k} \right) + \delta_{0} \right) \right]^{2} $$ を最小化する最適解 $l_{1} , \cdots , l_{n} , \delta_{0}$ を見つけることになる。

固有的定常性の定義: ユークリッド空間の固定された部分集合 $D \subset \mathbb{R}^{r}$ で、確率変数 $Y(s) : \Omega \to \mathbb{R}^{1}$ の集合である空間過程 $\left\{ Y(s) \right\}_{s \in D}$ と方向ベクトル $\mathbf{h} \in \mathbb{R}^{r}$ を考える。具体的に $n \in \mathbb{N}$ 個のサイト^siteを $\left\{ s_{1} , \cdots , s_{n} \right\} \subset D$ のように表し、$Y(s)$ はすべての $s \in D$ に対して分散が存在すると仮定する。$\left[ Y \left( s + \mathbf{h} \right) - Y(s) \right]$ の平均が $0$ でありながら分散が唯一 $\mathbf{h}$ にのみ依存する場合、$\left\{ Y \left( s_{k} \right) \right\}$ が固有的定常性^{intrinsic Stationarity}を持つと言われる。 $$ \begin{align*} E \left[ Y \left( s + \mathbf{h} \right) - Y(s) \right] =& 0 \\ \text{Var} \left[ Y \left( s + \mathbf{h} \right) - Y(s) \right] =& 2 \gamma ( \mathbf{h} ) \end{align*} $$

ここで $\left\{ Y \left( s_{k} \right) \right\}_{k=1}^{n}$ が固有的定常性を持つ場合、$\sum_{k} l_{k} = 1$ という制約条件を置くことにより、 $$ \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*} $$ となるようにできる。これにより、我々が最小化するべき目的関数は $\delta_{0}$ が外れた次の形になる。 $$ E \left[ Y \left( s_{0} \right) - \sum_{k} l_{k} Y \left( s_{k} \right) \right]^{2} + \delta_{0}^{2} $$ ここで $\delta_{0}^{2}$ は予測と関係ない。実際にもモデルが $\mathbf{Y} = \mu \mathbf{1} + \varepsilon$ であれば、$\delta_{0}$ は $\mu$ に該当し、$\delta_{0} = 0$ として $$ E \left[ Y \left( s_{0} \right) - \sum_{k} l_{k} Y \left( s_{k} \right) \right]^{2} $$ としても問題ない。今、$a_{0} = 1$ とし、$a_{k} = - l_{k}$ とすると、 $$ 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} $$ となるため、我々は次の最適化問題をラグランジュ乗数法で解くことになる。 $$ \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} $$

---

パート2. セミバリオグラム $\gamma$

今、$E \left[ \sum_{k=0}^{n} a_{k} Y \left( s_{k} \right) \right]^{2}$ を「我々がデータから計算できると仮定できる」セミバリオグラムに依存した形で表してみよう。

セミバリオグラムの定義: ユークリッド空間の固定された部分集合 $D \subset \mathbb{R}^{r}$ で、確率変数 $Y(s) : \Omega \to \mathbb{R}^{1}$ の集合である空間過程 $\left\{ Y(s) \right\}_{s \in D}$ と方向ベクトル $\mathbf{h} \in \mathbb{R}^{r}$ を考える。具体的に $n \in \mathbb{N}$ 個のサイトを $\left\{ s_{1} , \cdots , s_{n} \right\} \subset D$ のように表し、$Y(s)$ はすべての $s \in D$ に対して分散が存在すると仮定する。次のように定義される $2 \gamma ( \mathbf{h} )$ をバリオグラム^variogramと呼ぶ。 $$ 2 \gamma ( \mathbf{h} ) := E \left[ Y \left( s + \mathbf{h} \right) - Y(s) \right]^{2} $$ 特にバリオグラムの半分 $\gamma ( \mathbf{h} )$ をセミバリオグラム^{semivariogram}と呼ぶ。

$\gamma \left( s_{i} - s_{j} \right)$ を2つのサイト $s_{i}, s_{j}$ の間の方向ベクトルによるセミバリオグラムとしよう。$\sum_{0=1}^{n} a_{k} = 0$ を満たす任意の集合 $\left\{ a_{k} : k = 1 , \cdots , n \right\} \subset \mathbb{R}$ に対して次が成立する。 $$ \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} $$ これは次の展開 $$ \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} \text{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*} $$ で確認できる。今、$\gamma_{ij} = \gamma \left( s_{i} - s_{j} \right)$ とし、$\gamma_{0j} = \gamma \left( s_{0} - s_{j} \right)$ とすれば、我々の目的関数は $$ \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*} $$ のように表され、ラグランジュ乗数法によって制約条件 $\sum_{i} l_{i} = 1$ にラグランジュ乗数^{lagrange Multiplier}を掛けて引くことで、 $$ \ - \sum_{i} \sum_{j} l_{i} l_{j} \gamma_{ij} + 2 \sum_{i} l_{i} \gamma_{0i} - \lambda \sum_{i} l_{i} $$ を得る。したがって、$l_{i}$ に対して偏微分して $$ \ - \sum_{j} l_{j} \gamma_{ij} + \gamma_{0i} - \lambda = 0 $$ となる時、目的関数を最小化できることがわかる。

---

パート3. 最適解

ここで、具体的な最適解の公式を導出する。行列 $\Gamma \in \mathbb{R}^{n \times n}$ の $(i,j)$ 成分を ▷

eq74◁ とし、つまり $\left( \Gamma \right)_{ij} := \gamma_{ij}$ とし、ベクトル $\gamma_{0} \in \mathbb{R}^{n}$ を $\left( \gamma_{0} \right)_{i} := \gamma_{0i}$ のように定義しよう。係数のベクトルも同様に $l := \left( l_{1} , \cdots , l_{n} \right) \in \mathbb{R}^{n}$ とすると、パート2で得た式は次のような行列/ベクトル形式で表すことができる。 $$ \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*} $$ 一方、制約条件で $$ \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 $$ となるので、$\mathbf{1}^{T} l = 1$ も得られる。ここで、$\mathbf{x}^{T}$ は $\mathbf{x}$ の転置を表す。まず単独で $\lambda$ は $$ \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*} $$ となり、整理すると $$ \ - \lambda = {{ 1 - \mathbf{1}^T \Gamma^{-1} \gamma_{0} } \over { \mathbf{1}^T \Gamma^{-1} \mathbf{1} }} $$ のように表すことができる。今、$l$ はほぼ求めたも同然である。 $$ \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*} $$

■