Go to Qiita Advent Calendar 2024 Top
LoginSignup
7
7

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@faux(Yuma Uchiumi)

Sequential Gaussian Simulation: 少ない観測データから空間ヒートマップを作成する

Last updated at Posted at 2020-06-08

はじめに

  • このポストでは、空間統計モデルSequetial Gaussian Simulationについて解説します。
  • Googleで検索した限り、日本語の文献がほとんどなかったので、自分の学習記録も兼ねてここに記します。
  • まず、空間統計に固有の概念である「クリギング」「バリオグラム」について述べ、次に、実用的な空間回帰の手法として「Sequetial Gaussian Simulation」について説明します。

Cf.) Sequetial Gaussian Simulationの元論文
"Geostatistics in soil science: state-of-the-art and perspectives" (P.Goovaerts, 1999)

空間統計学(Geostatistics)

空間統計とは？

※ 空間統計学への招待から引用。

空間統計(Geostatistics)とは、ベクトル空間（$\mathbb{R}^2$,$\mathbb{R}^3$）中に生起するランダムな現象の統計解析の総称。

空間統計で扱われるモデルは以下の3つに大別される。らしい。
※参考：空間統計学への招待

  • 点パターン（point pattern）
    • 空間のランダムな点の集まり（点の位置と個数が共にランダム）のモデル。生態学（ecology）、森林学（forestry）、天文学、地理学（geography）等で使われる。
    • （例） 森林地帯の木の位置、星の位置、ある州の町の位置、町の中の八百屋の位置
  • ランダム閉集合（RACS、 random closed set）
    • ある空間の閉集合の全体が作る空間上の確率分布。ステレオロジー（stereology）やモルフォロジー（morphology）といった手法を用いて解析する。前者は、組織の平均特性量をその断面の観測値から推定する方法で、解剖学や鉱山学で用いられている。後者は、鉱石や生物組織断面にあらわれる複雑なパターンを、集合論的演算で変形（平滑化）することで特性量を抜き出す方法である。
  • 空間回帰（spatial regression）
    • 確率場（時系列の多次元版）に近い。観測値の分布には正規分布を仮定する。クリギングなどがある。
    • クリギング: 南アフリカの鉱山技術者Krigeが考案した方法。複数個のボーリングのデータ（サンプル）から採石場全体の鉱山の総含有量を推定する。発祥は鉱山学だが、気象学、森林学、農学、地図作成、気候（風土）学、水産学などに応用されている。

空間統計の歴史

空間統計学(Geostatistics)は、測定技術(GPS)や情報技術(BigData)の発展によってはじまった統計学の新しい領域といえます。観測される地理データそのものが時空間的であり、物理現象を直接とらえていることが魅力的。
以下、空間統計学への招待から引用。

空間統計は、いろいろな学問で必要に応じて開発されたデータ解析手法の総称である。これらの手法は古くから用いられてきたが、それらが一つの理論体系として記述されたのは1970年代から1980年代の間である。そして近年、空間統計が注目を集めるようになった背景には、計算機や計測機器の発達がある。空間統計で扱うような大規模なデータが蓄積され、また解析できるようになった為である。
（例）人口衛星による観測データ（画像）、CTスキャンによる画像（場合によっては3D）

データの空間的相関(Spatial Correlation)

バリオグラム(Matheron,1963)は、ある観測変数$Z$の"空間的な自己相関"つまり「データが距離と方向にどのような関係を持つか」を測る量です。

※似たような量として、"時間的な自己相関"を測る「コレログラム(Correlogram)」"空間的な相互相関"を測る「コバリオグラム(Covariogram)」などがあります。

自己相関(auto-correlation) 相互相関(cross-correlation)
空間の変化 バリオグラム(Variogram) コバリオグラム(Covariogram)
時間の変化 コレログラム(Correlogram) ココレログラム(Cocorrelogram)

バリオグラム(Variogram)の定義

※参考：「バリオグラム」（小林景, 慶應義塾大学）

平面上の各点${\bf x} \in \mathbb{R}^2$に，確率的な値$Z({\bf x})$が割り当てられているとします（こういったものを確率場という）。気象データの例でいうと、$\bf x$は必ずしも観測地点とは限らない日本地図上の各点であり，$Z({\bf x})$はそこでの降雨量などの気象値です。つまり、観測地点以外での気象値も含めて統計的推測を行います。このとき、バリオグラム$\gamma (\geq 0)$は2地点${\bf x}, {\bf x}'$における$Z$の観測値の差の分散として定義されます。

$$
\gamma({\bf x}, {\bf x}') := Var(Z({\bf x}) - Z({\bf x}')), ~~~ \forall {\bf x}, {\bf x}' \in \mathbb{R}^2
$$

バリオグラム$\gamma$は、$Z$の(空間的)共分散を用いて以下のように書き換えることもできます。

$$
\gamma({\bf x}, {\bf x}') = Var(Z({\bf x})) + Var(Z({\bf x}')) - 2Cov(Z({\bf x}), Z({\bf x}'))
$$

バリオグラム(Variogram)の性質

観測データZの定常性

観測変数$Z$に対して

定常性:「任意の${\bf x},{\bf x}' \in \mathbb{R}^2$に対して、$Z({\bf x})$と$Z({\bf x}')$は同じ確率分布に従う」

を仮定すると、2地点${\bf x}, {\bf x}'$における$Y$の相関係数

$$
Corr(Z({\bf x}), Z({\bf x}'))
:= \frac{Cov(Z({\bf x}), Z({\bf x}'))}{\sqrt{Var(Z({\bf x}))Var(Z({\bf x}'))}}
$$

が各点での分散$Var(Z({\bf x}))$は$\bf x$に依存しなくなるため，これを単に$Var(Z)$と表すと以下のようにバリオグラムと相関係数の関係が求まります。

（定常性の仮定下での）バリオグラムと空間的相関との関係:

$$
\begin{align}
Corr(Z({\bf x}), Z({\bf x}'))
&= \frac{Cov(Z({\bf x}), Z({\bf x}'))}{Var(Z)}
= 1 - \frac{\gamma({\bf x}, {\bf x}')}{2 Var(Z)}
\end{align}
$$

$$
\begin{align}
\gamma({\bf x}, {\bf x}')
&= 2 Var(Z) - 2Cov(Z({\bf x}), Z({\bf x}')) \
&= 2 Var(Z) - 2Var(Z) Corr(Z({\bf x}), Z({\bf x}')) \
&= 2 Var(Z) (1 - Corr(Z({\bf x}), Z({\bf x}')) )
\end{align}
$$

つまり、定常性を仮定すると、2地点間の空間的相関 $Cov(Z({\bf x}), Z({\bf x}'))$ が $1$ に近づくとバリオグラム $\gamma({\bf x}, {\bf x}')$は$0$へ近づき、逆に2地点間の空間的相関 $Cov(Z({\bf x}), Z({\bf x}'))$ が $0$ に近づくと一定値の空間分散 $2Var(Z)$ に近づきます。

観測データZの等方性

観測変数$Z$に対して

等方性(isotropy): 「任意の${\bf x}, \bf{h} \in \mathbb{R}^2$に対して、$Z({\bf x} + {\bf h}) - Z({\bf x})$の確率分布が$r = || {\bf h} ||$（ノルム,距離）にのみ依存する」

を仮定すると、2地点${\bf x}, {\bf x}'$における$Y$のバリオグラム$\gamma({\bf x},{\bf x}')$は、$r = ||{\bf x}' - {\bf x}||$を用いて、$\gamma(r)$のように書けます。

※ 等方性(isotropy)を持たないことを異方性(anisotropy)という

データの空間的補間（Spatial Interpolation）

GPSなどの計測器から実際に得られる地理情報・空間データでは、平面上の任意の地点${\bf x} \in \mathbb{R}^2$における観測値$Z({\bf x})$が得られるわけではありません。そのため、平面上にある$N$コの観測地点$({\bf x}_1, \dots, {\bf x}_N)$で得られた$N$コの観測データ$(Z({\bf x}_1), \dots, Z({\bf x}_N))$を用いて

  1. 任意の2地点${\bf x}, {\bf x}'$における$Z$のバリオグラム$\gamma({\bf x},{\bf x}')$
  2. 任意の1地点${\bf x}$における$Z$の値$Z({\bf x})$

を推定する問題を考えます。バリオグラム$\gamma({\bf x},{\bf x}')$が求められれば、$Z$の空間分布(=ヒートマップ)が得られます。ヒートマップから任意の${\bf x}$における$Z$の値$Z({\bf x})$そのものを直接予測することもできます。

  • 空間的補間 or 内挿（spatial inteτpolation）：観測を行った地域内でデータのない地点の値を推定
  • 空間的補外 or 外挿（spatial extrapolation）：観測を行った地域外でデータのない地点の値を推定

近隣地点のデータのみを使った当てはめ法（局地的補間法）としては、スプライン平滑化・移動平均・クリギングが使われます。

クリギング(Kringing)

クリギング(Kringing)とは、空間上の限られた地点での観測値から、任意の地点での値を予測する（＝空間的補間）ための手法である。

クリギング(Kringing): 南アフリカの鉱山技術者Krigeが考案した方法。複数個のボーリングのデータ（サンプル）から採石場全体の鉱山の総含有量を推定する。発祥は鉱山学だが、気象学、森林学、農学、地図作成、気候（風土）学、水産学などに応用されている。（※ 空間統計学への招待から引用）

クリギングでは、定常性を仮定した次のような条件が成立することが仮定されます。すなわち, 任意の2地点間の$Z$の差(確率変数)の期待値は$0$で、分散は地点間の距離$||{\bf x}' - {\bf x}||$のみに依存すると仮定されます。

$$
\begin{align}
\mathbb{E}[Z({\bf x}) - Z({\bf x}')] &= 0 \
Var(Z({\bf x}) - Z({\bf x}')) &= \gamma({\bf x}, {\bf x}')
\end{align}
$$

$\gamma$を推定し、それを用いて空間補間することで$Z$の予測を行うのがクリギングである。

普通クリギング（Ordinary Kriging）

クリギングでは、$Z({\bf x}_0)$ に対する推定値を求めるために、その周囲の$N$個の観測値 $\{ Z({\bf x}_1), \dots, Z({\bf x}_N) \}$ の重み付き平均を使います。すなわち、ある地点 ${\bf x}_0$ の値 $Z({\bf x}_0)$ に対する推定値 $\hat{Z}({\bf x}_0)$ は

$$
Z({\bf x}_0) \approx \hat{Z}({\bf x}_0) = \sum_{i=1}^{N} w_i Z({\bf x}_i), ~~~
s.t. ~~~ \sum_{i=1}^{N} w_i = 1
$$

となります。クリギングが通常の重み付き平均法と異なる点は、重み（パラメータ）$\{ w_1, \dots, w_N \}$ を求める際に、誤差ではなく観測値から得られたバリオグラムを利用することです。（真値 $Z({\bf x}_0)$ がわからないため、推定誤差 $\hat{Z}({\bf x}_0) - Z({\bf x}_0)$ は計算できないため。）

以下では、最も広く利用されている普通クリギング（Ordinary Kriging）を紹介します。普通クリギングでは、理論バリオグラム $\gamma(h)$ の代わりに、標本バリオグラム $\gamma^{(OK)}(h)$ を使います。

理論バリオグラム
$$
\gamma(h) = \frac{1}{2} ~ \mathbb{E} \left[ { \left\{ Z({\bf x} + {\bf h}) - Z({\bf x}) \right\} }^2 \right]
$$
標本バリオグラム（経験バリオグラム）
$$
\begin{align}
\gamma^{(OK)}(h)
&= \frac{1}{2 |N(h)|} \sum_{(i,j) \in N(h)} {\left( Z({\bf x}_i) - Z({\bf x}_j) \right) }^2 \
&= \frac{1}{2N} \sum { \left\{ Z({\bf x}_i + {\bf h}) - Z({\bf x}_i) \right\} }^2 \
\end{align}
$$

推定誤差に対する平均・分散を計算すると、次のようになります。

$$
\begin{align}
\text{(}Z_0\text{の推定誤差の平均)}
&= \mathbb{E} [ \hat{Z}({\bf x}_0) - Z({\bf x}_0) ] \
&= \mathbb{E} \left[ \sum_{i=1}^{N} w_i Z({\bf x}_i) - Z({\bf x}_0) \right] \
&= \mathbb{E} \left[ \sum_{i=1}^{N} w_i Z({\bf x}_i) - Z({\bf x}_0) \left( \sum_{i=1}^{N} w_i \right) \right] \
&= \sum_{i=1}^{N} w_i ~ \mathbb{E}\left[ Z({\bf x}_i) - Z({\bf x}_0) \right] \
&= 0
\\ \
\text{(}Z_0\text{の推定誤差の分散)}
&= \mathbb{E} \left[ { \left\{ \hat{Z}({\bf x}_0) - Z({\bf x}_0) \right\} }^2 \right] \
&= - \gamma({\bf x}_0, {\bf x}_0) - \sum_{i=1}^{N} \sum_{j=1}^{N} w_i w_j \gamma({\bf x}_i, {\bf x}_j) + 2 \sum_{i=1}^{N} w_i \gamma({\bf x}_i, {\bf x}_0) \
&= - \gamma^{(OK)}(0) - \sum_{i=1}^{N} \sum_{j=1}^{N} w_i w_j \gamma^{(OK)}( || {\bf x}_i - {\bf x}_j || ) + 2 \sum_{i=1}^{N} w_i \gamma^{(OK)}(|| {\bf x}_i - {\bf x}_0 ||) \
&= \sigma^2_{E}( \hat{Z}({\bf x}_0) ~\vert~ w_1, \dots, w_N )
\end{align}
$$

よって、“最良”の重み $\{ w_1, \dots, w_N \}$ を求めるために，推定誤差の分散$\sigma^2_{E}$をできる限り小さくする最適化問題を解く必要があります。

$$
\begin{align}
&\underset{ \{ w_1, \dots, w_N \} }{\min}
~ \sigma^2_{E}( \hat{Z}({\bf x}_0) ~\vert~ w_1, \dots, w_N ) \
& ~~~~ s.t. ~~~ \sum_{i=1}^{N} w_i = 1
\end{align}
$$

シミュレーション

ここまで、データの空間的相関としての「バリオグラム」と、バリオグラムを評価指標として用いた空間回帰「クリギング」について考えてきました。ここでは代表的なシミュレーション手法であるSequential Gaussian Simiulation (SGSIM, SGS)について説明します。

Sequential Gaussian Simulation

SGSIMは、普通クリギングとモンテカルロ法の組み合わせで実行できるシュミレーション手法で、最も一般的な手法です。すなわち、SGSIMでは

  1. 観測データ$\{ Z({\bf x}_1), \dots, Z({\bf x}_N) \}$に対する普通クリギングによる$\gamma$と空間分布の推定
  2. Monte Carloサンプリングによる尤もらしい疑似データ点$Z({\bf x}_0)$の追加

を逐次的(Sequential)に繰り返します。SGSIMの”sequential”は時系列性ではなく、上記の1と2を逐次的に繰り返すという意味です。つまり、SGSIMでは時間情報を扱えっておらず。純粋に空間的相関を推定しています。SGSIMはある決められた時刻での「ある変数(人口、感染率、降雨量)のヒートマップを計算する」目的に対して有用で、広く使われているモデルです。

  • SGSIMの長所
    • 比較的シンプルなモデル、実装が容易
    • 少ない観測点からでも、MCサンプリングにより空間上の任意点における推定値を計算できる
  • SGSIMの短所
    • MCサンプリングの計算量
    • 時間情報を考慮していない。（例: トレンド、季節変動、周期性など）SGSIMの発展モデルもいろいろある。

SGSIMの逐次更新アルゴリズム

SGSIMの逐次更新アルゴリズムについて述べます。

Step 0 (正規化) 観測データ$ \left[ Z({\bf x}_1),\dots,Z({\bf x}_N) \right]$を標準化する。すなわち、標本平均が0、標本分散が1になるように値のスケールを調節する。（z-score transformation）

$$
\begin{align}
&Z({\bf x}_i) \gets \frac{Z({\bf x}_i) - \mu_{Z} }{\sigma^2_{Z}} \
where ~~~ \mu_Z =& \frac{1}{N} \sum_{j=1}^{N} Z({\bf x}_j)
, ~~~
\sigma^2_{Z} = \frac{1}{N} \sum_{j=1}^{N} { \{ Z({\bf x}_j) - \mu_Z \} }^2
\end{align}
$$

Step 1 (Kringing) Kringing modelを使って空間回帰を行い、推定値(平均)と分散を求める。

$$
\begin{align}
\hat{Z}({\bf x})
&= \sum_{i=1}^{N} \beta_i Z({\bf x}_i) \
\sigma^2({\bf x})_{SK}
&= C({\bf 0}) - \sum_{i=1}^{N} \alpha_i Cov({\bf x}_i, {\bf x})
\end{align}
$$

Step 2 (Resampling) 正規乱数 $R({\bf x}) \sim \mathcal{N}({\bf 0}, \sigma^2({\bf x})_{SK})$ を加える
$$
\hat{Z}_s({\bf x}) = \hat{Z}({\bf x}) + R({\bf x}) \
$$
つまり、$\hat{Z}_s({\bf x})$は、$Z({\bf x})$に対する予測分布$\hat{Z}_s({\bf x}) \sim \mathcal{N}(\hat{Z}({\bf x}), \sigma^2({\bf x})_{SK})$からサンプリングされた値。
4. 観測データ$ \left[ Z({\bf x}_1),\dots,Z({\bf x}_N), \hat{Z}_s({\bf x}) \right]$を標準化する。

クリギングをして平均(kringing mean)と分散(kringing variance)を求める

参考文献

7
7
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
7
7

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?