ユニバーサル・クリギング
📂統計的分析 ユニバーサル・クリギング モデル 空間データ分析 で、ランダムフィールド Y = ( Y ( s 1 ) , ⋯ , Y ( s n ) ) \mathbf{Y} = \left( Y \left( s_{1} \right) , \cdots , Y \left( s_{n} \right) \right) Y = ( Y ( s 1 ) , ⋯ , Y ( s n ) ) の平均 μ ∈ R \mu \in \mathbb{R} μ ∈ R と共分散行列 Σ ∈ R n × n \Sigma \in \mathbb{R}^{n \times n} Σ ∈ R n × n と多変量正規分布 ε ∼ N n ( 0 , Σ ) \varepsilon \sim N_{n} \left( \mathbf{0} , \Sigma \right) ε ∼ N n ( 0 , Σ ) に従うとき、次のモデル
Y = μ 1 + ε
\mathbf{Y} = \mu \mathbf{1} + \varepsilon
Y = μ 1 + ε
を使って新しいサイトsite s 0 s_{0} s 0 のY ( s 0 ) Y \left( s_{0} \right) Y ( s 0 ) の推定値をオーディナリークリギング推定値 ordinary Kriging Estimate と呼ぶ。このようにしてモデルを設定し推定する行為自体もクリギング と呼ばれる。
ユニバーサルクリギング X : = [ 1 X 1 ( s 1 ) ⋯ X p ( s 1 ) 1 X 1 ( s 2 ) ⋯ X p ( s 2 ) ⋮ ⋮ ⋱ ⋮ 1 X 1 ( s n ) ⋯ X p ( s n ) ]
\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 1 ⋮ 1 X 1 ( s 1 ) X 1 ( s 2 ) ⋮ X 1 ( s n ) ⋯ ⋯ ⋱ ⋯ X p ( s 1 ) X p ( s 2 ) ⋮ X p ( s n )
独立変数 X 1 , ⋯ , X p X_{1} , \cdots , X_{p} X 1 , ⋯ , X p に対して、デザインマトリックス と回帰係数のベクトル β = ( β 0 , β 1 , ⋯ , β p ) \beta = \left( \beta_{0} , \beta_{1} , \cdots , \beta_{p} \right) β = ( β 0 , β 1 , ⋯ , β p ) を使った多変量正規分布 ε ∼ N n ( 0 , Σ ) \varepsilon \sim N_{n} \left( \mathbf{0} , \Sigma \right) ε ∼ N n ( 0 , Σ ) に従う次のモデル
Y = X β + ε
\mathbf{Y} = \mathbf{X} \beta + \varepsilon
Y = X β + ε
を通じ、新しいサイトsite s 0 s_{0} s 0 のx ( s 0 ) = ( x 1 ( s 0 ) , ⋯ , x p ( s 0 ) ) \mathbf{x} \left( s_{0} \right) = \left( x_{1} \left( s_{0} \right) , \cdots , x_{p} \left( s_{0} \right) \right) x ( s 0 ) = ( x 1 ( s 0 ) , ⋯ , x p ( s 0 ) ) を使用しY ( s 0 ) Y \left( s_{0} \right) Y ( s 0 ) の推定を行うことをユニバーサルクリギング という。
1 = ( 1 , ⋯ , 1 ) \mathbf{1} = (1 , \cdots , 1) 1 = ( 1 , ⋯ , 1 ) は、全ての成分が1 1 1 である1ベクトル 。N n ( 0 , Σ ) N_{n} \left( \mathbf{0} , \Sigma \right) N n ( 0 , Σ ) は多変量正規分布 を意味する。説明 ユニバーサルクリギング とは、簡単に言えば、クリギング に重回帰分析 が加えられたクリギングのことで、オーディナリークリギングのμ \mu μ はユニバーサルクリギングのβ 0 \beta_{0} β 0 と同じ役割を果たす。時系列分析 としてみれば、主に興味を持っている変数Y Y Y を説明するための情報として、純粋に空間的な構造だけでなく、他の変数も考慮に入れられる動的回帰モデル とみなせる。
式を見ると、実際には回帰分析の部分が異なって見えるかもしれないが、X β \mathbf{X} \beta X β を左辺に移動して
( Y − X β ) = 0 + ε
\left( \mathbf{Y} - \mathbf{X} \beta \right) = \mathbf{0} + \varepsilon
( Y − X β ) = 0 + ε
とした場合、それは単に平均ベクトルがゼロベクトルμ = 0 \mu = \mathbf{0} μ = 0 のオーディナリークリギング そのものだ。
公式 ランダムフィールド { Y ( s k ) } k = 1 n \left\{ Y (s_{k}) \right\}_{k=1}^{n} { Y ( s k ) } k = 1 n が固有の定常空間プロセス だとし、新しく予測したい地点をs 0 s_{0} s 0 とする。バリオグラム 2 γ 2 \gamma 2 γ に関してベクトル γ ∈ R n \gamma \in \mathbb{R}^{n} γ ∈ R n を次のように定義する。あるベクトルλ = ( λ 1 , ⋯ , λ n ) \lambda = \left(\lambda_{1} , \cdots ,\lambda_{n} \right) λ = ( λ 1 , ⋯ , λ n ) についてY ( s 0 ) Y \left( s_{0} \right) Y ( s 0 ) の最良線形不偏予測量bLUP, Best Linear Unbiased Predictor は
λ T Y = [ λ 1 ⋯ λ n ] [ Y ( s 1 ) Y ( s 2 ) ⋮ Y ( s n ) ] = ∑ k = 1 n λ k Y ( s k )
\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)
λ T Y = [ λ 1 ⋯ λ n ] Y ( s 1 ) Y ( s 2 ) ⋮ Y ( s n ) = k = 1 ∑ n λ k Y ( s k )
であり、ベクトルλ \lambda λ は具体的に次のように求められる。
λ = Σ − 1 γ + Σ − 1 X ( X T Σ − 1 X ) − 1 ( x 0 − X T Σ − 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 X ( X T Σ − 1 X ) − 1 ( x 0 − X T Σ − 1 γ )
導出 パート1. 条件付き期待値
E [ ( Y ( s 0 ) − h ( y ) ) 2 ∣ y ]
E \left[ \left. \left( Y \left( s_{0} \right) - h \left( \mathbf{y} \right) \right)^{2} \right| \mathbf{y} \right]
E [ ( Y ( s 0 ) − h ( y ) ) 2 y ]
目標は、データy = ( y ( s 1 ) , ⋯ , y ( s n ) ) \mathbf{y} = \left( y \left( s_{1} \right) , \cdots , y \left( s_{n} \right) \right) y = ( y ( s 1 ) , ⋯ , y ( s n ) ) が与えられた時に上記の誤差を最小化するh ( y ) h \left( \mathbf{y} \right) h ( y ) を求めることである。
E [ ( Y ( s 0 ) − h ( y ) ) 2 ∣ y ] = E [ ( Y ( s 0 ) − E ( Y ( s 0 ) ∣ y ) ) 2 ∣ y ] + [ E ( Y ( s 0 ) ∣ y ) − h ( y ) ] 2 ≥ E [ ( Y ( s 0 ) − E ( Y ( s 0 ) ∣ 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*}
= ≥ E [ ( Y ( s 0 ) − h ( y ) ) 2 y ] E [ ( Y ( s 0 ) − E ( Y ( s 0 ) ∣ y ) ) 2 y ] + [ E ( Y ( s 0 ) ∣ y ) − h ( y ) ] 2 E [ ( Y ( s 0 ) − E ( Y ( s 0 ) ∣ y ) ) 2 y ]
この式が成立するための必要十分条件はh ( y ) = E ( Y ( s 0 ) ∣ y ) h \left( \mathbf{y} \right) = E \left( Y \left( s_{0} \right) | \mathbf{y} \right) h ( y ) = E ( Y ( s 0 ) ∣ y ) である。
パート2. 条件付き多変量正規分布
基本的に、与えられたデータ は多変量正規分布 N n ( 0 , Σ ) N_{n} \left( \mathbf{0} , \Sigma \right) N n ( 0 , Σ ) に従うという前提の下、新たなY 0 : = Y ( S 0 ) ∼ N ( μ 0 , σ 0 2 ) Y_{0} := Y \left( S_{0} \right) \sim N \left( \mu_{0} , \sigma_{0}^{2} \right) Y 0 := Y ( S 0 ) ∼ N ( μ 0 , σ 0 2 ) を追加しても依然として( 1 + n ) (1+n) ( 1 + n ) 変量の正規分布に従うという仮定から始める。これをブロック行列 形式で表すと次のようになる。
[ Y 0 Y 1 ⋮ Y n ] = [ Y 0 Y ] ∼ N 1 + n ( [ μ 0 μ ] , [ σ 0 2 γ 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)
Y 0 Y 1 ⋮ Y n = [ Y 0 Y ] ∼ N 1 + n ( [ μ 0 μ ] , [ σ 0 2 γ γ T Σ ] )
多変量正規分布の条件付き平均と分散 :
X = [ X 1 X 2 ] : Ω → R n μ = [ μ 1 μ 2 ] ∈ R n Σ = [ Σ 11 Σ 12 Σ 21 Σ 22 ] ∈ R n × 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 = μ = Σ = [ X 1 X 2 ] [ μ 1 μ 2 ] [ Σ 11 Σ 21 Σ 12 Σ 22 ] : Ω → R n ∈ R n ∈ R n × n
ジョルダンブロック形式 で表されたX \mathbf{X} X 、μ \mu μ 、Σ \Sigma Σ に対し多変量正規分布 に従うランダムベクトル X ∼ N n ( μ , Σ ) \mathbf{X} \sim N_{n} \left( \mu , \Sigma \right) X ∼ N n ( μ , Σ ) が与えられているとする。すると、条件付き確率ベクトル X 1 ∣ X 2 : Ω → R m \mathbf{X}_{1} | \mathbf{X}_{2} : \Omega \to \mathbb{R}^{m} X 1 ∣ X 2 : Ω → R m も依然として多変量正規分布に従い、具体的には次のように母平均ベクトルと母共分散行列 を持つ。
X 1 ∣ X 2 ∼ N m ( μ 1 + Σ 12 Σ 22 − 1 ( X 2 − μ 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)
X 1 ∣ X 2 ∼ N m ( μ 1 + Σ 12 Σ 22 − 1 ( X 2 − μ 2 ) , Σ 11 − Σ 12 Σ 22 − 1 Σ 21 )
したがって、データY \mathbf{Y} Y が与えられた時のY 0 Y_{0} Y 0 の期待値 はs 0 s_{0} s 0 でのx 0 = x ( s 0 ) ∈ R p \mathbf{x}_{0} = \mathbf{x} \left( s_{0} \right) \in \mathbb{R}^{p} x 0 = x ( s 0 ) ∈ R p に対して
E ( Y ( s 0 ) ∣ Y ) = μ 0 + γ T Σ − 1 ( Y − X β ) = x 0 T β + γ 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*}
= = E ( Y ( s 0 ) ∣ Y ) μ 0 + γ T Σ − 1 ( Y − X β ) x 0 T β + γ T Σ − 1 ( Y − X β )
である。一方で、回帰係数のベクトル β \beta β は、多変量正規分布に従うε \varepsilon ε の各成分が独立でないため、単純なOLSordinary Least Square のβ ^ = ( X T X ) − 1 X T Y \hat{\beta} = \left( \mathbf{X}^{T} \mathbf{X} \right)^{-1} \mathbf{X}^{T} \mathbf{Y} β ^ = ( X T X ) − 1 X T Y ではなく、WLSweighted Least Square のβ ^ = ( X T Σ − 1 X ) − 1 X T Σ − 1 Y \hat{\beta} = \left( \mathbf{X}^{T} \Sigma^{-1} \mathbf{X} \right)^{-1} \mathbf{X}^{T} \Sigma^{-1} \mathbf{Y} β ^ = ( X T Σ − 1 X ) − 1 X T Σ − 1 Y が使用される必要がある。
E ( Y ( s 0 ) ∣ Y ) = x 0 T β + γ T Σ − 1 Y − γ T Σ − 1 X β = x 0 T β + γ T Σ − 1 Y − γ T Σ − 1 X β = γ T Σ − 1 Y + x 0 T ( X T Σ − 1 X ) − 1 X T Σ − 1 Y − γ T Σ − 1 X ( X T Σ − 1 X ) − 1 X T Σ − 1 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*}
= = = E ( Y ( s 0 ) ∣ Y ) x 0 T β + γ T Σ − 1 Y − γ T Σ − 1 X β x 0 T β + γ T Σ − 1 Y − γ T Σ − 1 X β γ T Σ − 1 Y + x 0 T ( X T Σ − 1 X ) − 1 X T Σ − 1 Y − γ T Σ − 1 X ( X T Σ − 1 X ) − 1 X T Σ − 1 Y
ここでΣ \Sigma Σ は対称行列 であるため( Σ − 1 ) T = Σ − 1 \left( \Sigma^{-1} \right)^{T} = \Sigma^{-1} ( Σ − 1 ) T = Σ − 1 であり、Y \mathbf{Y} Y と組み合わせると次のようなλ T Y \lambda^{T} \mathbf{Y} λ T Y 形式を得る。
γ T Σ − 1 Y + x 0 T ( X T Σ − 1 X ) − 1 X T Σ − 1 Y − γ T Σ − 1 X ( X T Σ − 1 X ) − 1 X T Σ − 1 Y = [ Σ − 1 γ + Σ − 1 X ( X T Σ − 1 X ) − 1 ( x 0 − X T Σ − 1 γ ) ] T Y = λ T 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*}
= = γ T Σ − 1 Y + x 0 T ( X T Σ − 1 X ) − 1 X T Σ − 1 Y − γ T Σ − 1 X ( X T Σ − 1 X ) − 1 X T Σ − 1 Y [ Σ − 1 γ + Σ − 1 X ( X T Σ − 1 X ) − 1 ( x 0 − X T Σ − 1 γ ) ] T Y λ T Y
■
実際の分析では、σ 0 2 \sigma_{0}^{2} σ 0 2 とΣ \Sigma Σ の推定値は、セミバリオグラムモデル を通じて次のように定められる。
σ 0 2 = σ 2 + τ 2 Σ = σ 2 H ( ϕ ) + τ 2 I
\begin{align*}
\sigma_{0}^{2} =& \sigma^{2} + \tau^{2}
\\ \Sigma =& \sigma^{2} H \left( \phi \right) + \tau^{2} I
\end{align*}
σ 0 2 = Σ = σ 2 + τ 2 σ 2 H ( ϕ ) + τ 2 I
ここで、τ 2 \tau^{2} τ 2 はナゲット効果の分散 (理論とは異なり、実際のデータ で距離に関係なく基本的に見られる共分散性)であり、I I I は単位行列 である。