📂統計的分析

ユニバーサル・クリギング

モデル

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

空間データ分析で、ランダムフィールド $\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{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} $$ 独立変数 $X_{1} , \cdots , X_{p}$ に対して、デザインマトリックスと回帰係数のベクトル$\beta = \left( \beta_{0} , \beta_{1} , \cdots , \beta_{p} \right)$を使った多変量正規分布$\varepsilon \sim N_{n} \left( \mathbf{0} , \Sigma \right)$に従う次のモデル $$ \mathbf{Y} = \mathbf{X} \beta + \varepsilon $$ を通じ、新しいサイト^site $s_{0}$ の$\mathbf{x} \left( s_{0} \right) = \left( x_{1} \left( s_{0} \right) , \cdots , x_{p} \left( s_{0} \right) \right)$ を使用し$Y \left( s_{0} \right)$ の推定を行うことをユニバーサルクリギングという。

---

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

説明

ユニバーサルクリギングとは、簡単に言えば、クリギングに重回帰分析が加えられたクリギングのことで、オーディナリークリギングの$\mu$はユニバーサルクリギングの$\beta_{0}$と同じ役割を果たす。時系列分析としてみれば、主に興味を持っている変数$Y$を説明するための情報として、純粋に空間的な構造だけでなく、他の変数も考慮に入れられる動的回帰モデルとみなせる。

式を見ると、実際には回帰分析の部分が異なって見えるかもしれないが、$\mathbf{X} \beta$を左辺に移動して $$ \left( \mathbf{Y} - \mathbf{X} \beta \right) = \mathbf{0} + \varepsilon $$ とした場合、それは単に平均ベクトルがゼロベクトル$\mu = \mathbf{0}$のオーディナリークリギングそのものだ。

公式

ランダムフィールド $\left\{ Y (s_{k}) \right\}_{k=1}^{n}$が固有の定常空間プロセスだとし、新しく予測したい地点を$s_{0}$とする。バリオグラム$2 \gamma$に関してベクトル$\gamma \in \mathbb{R}^{n}$を次のように定義する。あるベクトル$\lambda = \left(\lambda_{1} , \cdots ,\lambda_{n} \right)$について$Y \left( s_{0} \right)$の最良線形不偏予測量^{bLUP, Best Linear Unbiased Predictor}は $$ \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$は具体的に次のように求められる。 $$ \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. 条件付き期待値

$$ E \left[ \left. \left( Y \left( s_{0} \right) - h \left( \mathbf{y} \right) \right)^{2} \right| \mathbf{y} \right] $$ 目標は、データ$\mathbf{y} = \left( y \left( s_{1} \right) , \cdots , y \left( s_{n} \right) \right)$が与えられた時に上記の誤差を最小化する$h \left( \mathbf{y} \right)$を求めることである。 $$ \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 \left( \mathbf{y} \right) = E \left( Y \left( s_{0} \right) | \mathbf{y} \right)$である。

---

パート2. 条件付き多変量正規分布

基本的に、与えられたデータは多変量正規分布 $N_{n} \left( \mathbf{0} , \Sigma \right)$に従うという前提の下、新たな$Y_{0} := Y \left( S_{0} \right) \sim N \left( \mu_{0} , \sigma_{0}^{2} \right)$を追加しても依然として$(1+n)$変量の正規分布に従うという仮定から始める。これをブロック行列形式で表すと次のようになる。 $$ \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) $$

多変量正規分布の条件付き平均と分散: $$ \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*} $$ ジョルダンブロック形式で表された$\mathbf{X}$、$\mu$、$\Sigma$に対し多変量正規分布に従うランダムベクトル$\mathbf{X} \sim N_{n} \left( \mu , \Sigma \right)$が与えられているとする。すると、条件付き確率ベクトル$\mathbf{X}_{1} | \mathbf{X}_{2} : \Omega \to \mathbb{R}^{m}$も依然として多変量正規分布に従い、具体的には次のように母平均ベクトルと母共分散行列を持つ。 $$ \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) $$

したがって、データ$\mathbf{Y}$が与えられた時の$Y_{0}$の期待値は$s_{0}$での$\mathbf{x}_{0} = \mathbf{x} \left( s_{0} \right) \in \mathbb{R}^{p}$に対して $$ \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$の各成分が独立でないため、単純なOLS^{ordinary Least Square}の$\hat{\beta} = \left( \mathbf{X}^{T} \mathbf{X} \right)^{-1} \mathbf{X}^{T} \mathbf{Y}$ではなく、WLS^{weighted Least Square}の$\hat{\beta} = \left( \mathbf{X}^{T} \Sigma^{-1} \mathbf{X} \right)^{-1} \mathbf{X}^{T} \Sigma^{-1} \mathbf{Y}$が使用される必要がある。 $$ \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$は対称行列であるため$\left( \Sigma^{-1} \right)^{T} = \Sigma^{-1}$であり、$\mathbf{Y}$と組み合わせると次のような$\lambda^{T} \mathbf{Y}$形式を得る。 $$ \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*} $$

■

実際の分析では、$\sigma_{0}^{2}$と$\Sigma$の推定値は、セミバリオグラムモデルを通じて次のように定められる。 $$ \begin{align*} \sigma_{0}^{2} =& \sigma^{2} + \tau^{2} \\ \Sigma =& \sigma^{2} H \left( \phi \right) + \tau^{2} I \end{align*} $$ ここで、$\tau^{2}$はナゲット効果の分散(理論とは異なり、実際のデータで距離に関係なく基本的に見られる共分散性)であり、$I$は単位行列である。