空間データのためのブートストラップ法(解説編)

東京工業大学/株式会社Nospare の栗栖です．
この記事では空間データ・時空間データ分析に関する最近の研究成果「高次元空間データに対するガウス近似と空間ワイルドブートストラップ」 (Kurisu, Kato and Shao(2024)) について紹介します．今回は主に論文で考える空間データの設定や空間ワイルドブートストラップの実装手続きについて解説します．空間ワイルドブートストラップを利用した具体的な分析例は(データ分析編) で紹介したいと思います．

概要

様々な地点において観測された降水量や気温，大気中の化学物質の濃度などは空間データと呼ばれ，そのようなデータを扱う統計学は空間統計学と呼ばれます．特に近年ではリモートセンシング (遠隔探査) の発達により，各地点に置いて多くの変数 (例えばオゾン濃度，$\text{PM}_{2.5}$，$\text{PM}_{10}$，一酸化窒素・一酸化炭素濃度など) が観測できるようになってきました．また年間降水量などに代表される多くの空間データは複数の時点で空間データの観測が行われ，そのようなデータは時空間データと呼ばれます．この研究では空間データ・時空間データのための統計解析手法を新たに提案しました．

空間データ・時空間データ

まずこの研究で扱う空間データ・時空間データの構造について説明します．

空間データ

モデル : $\boldsymbol{Y} = \{\boldsymbol{Y}(\boldsymbol{s}) = (Y_{1}(\boldsymbol{s}),\cdots,Y_{p}(\boldsymbol{s}))' \in \mathbb{R}^{p}: \boldsymbol{s} \in R_{n}\}$．
即ち，このモデルでは観測領域 $R_{n}(\subset \mathbb{R}^{d})$ 内の点 $\boldsymbol{s}$ において $p$ 個の変数が観測される．空間データは特に $R_{n}\subset \mathbb{R}^{2}$ の場合に対応します．

データ : $\{\boldsymbol{Y}(\boldsymbol{s}_{i})\}_{i=1}^{n}$．
即ち，$R_{n}$ 内の地点 $\boldsymbol{s}_{1}, \cdots, \boldsymbol{s}_{n}$ でモデル $\boldsymbol{Y}$ を離散観測する状況を考えます．

(ポイント)

観測地点 $\boldsymbol{s}_{i}$ は観測領域 $R_{n}$ 内で不等間隔に散らばっていてもOK．
データ観測の枠組みについて，観測地点は規則的に配置された格子点上に分布する，という流儀もあります．ただし，各格子点に画素が与えられた画像データのような場合を除き，空間データは物理的制約により観測地点が不等間隔になる場合が多いため，この研究では格子点上でデータが観測される場合も特別な場合として扱える上記の設定で考えることにします．

時空間データ

モデル : $\boldsymbol{Z} = \{Z(t,\boldsymbol{s}) \in \mathbb{R} : t \in T_{n}, \boldsymbol{s} \in R_{n}\}$
このモデルは $R_{n}$ 内の点 $\boldsymbol{s}$ において時点 $t \in T_{n} (\subset \mathbb{R}\ \text{or}\ \mathbb{Z})$ で1つの値をもちます．

データ : $\{Z(t_{j},\boldsymbol{s}_{k})\}_{1 \leq j \leq p, 1 \leq k \leq n}$．
即ち，$R_{n}$ 内の地点 $\boldsymbol{s}_{i}$ で時系列 $(Z(t_{1},\boldsymbol{s}_{i}), \cdots, Z(t_{p},\boldsymbol{s}_{i}))'$ を離散観測する状況を考えます. また各地点で共通の観測時点を $t_{1}<\cdots <t_{p}$ としておきます．

(ポイント)

観測地点 $\boldsymbol{s}_{i}$ は観測領域 $R_{n}$ 内で不等間隔に散らばっていてもOK．
観測時点 $t_{j}$ も観測期間 $T_{n}$ 内で不等間隔に並んでいてもOK．
各地点での観測期間が長くとれる場合は $t_{p} \to \infty$ かつ $p = p_{n} \to \infty$ ($n \to \infty$)．

時空間データの例

下の図はアメリカの中西部の6つの州 (ノースダコタ，サウスダコタ，ネブラスカ，ミネソタ，ウィスコンシン，アイオワ) において 1919年から1994年まで毎年降水量の観測データが記録されている気象台の配置図で，それぞれの点が気象台の位置に対応しています (気象台の数は101地点)．この図を見ると，観測地点は観測領域内で比較的密になっている地域もあればそうでない地域もあり，不等間隔に散らばっていることがわかります．特にこの例では毎年各地点での年間降水量が観測されているので，年ごとに空間データが得られていることになります．さらに複数の時点で空間データが得られているのでデータ全体としては時空間データとみることができます．

より具体的にはこの例では

$R_{n} = (\text{6つの州の面積})$，
観測地点 $\boldsymbol{s}_{k}$ はそれぞれの気象台の位置，
観測時点 $t_{j}$ は1919年から1994年までの各年，
観測値 $Z(t_{j},\boldsymbol{s}_{k}) = (\text{地点 $\boldsymbol{s}_{k}$ の気象台における時点 $t_{j}$ での降水量の観測値})$

という対応になっています．各 $t_{j}$ を固定してみると $Y_{j}(\boldsymbol{s}_{i}) = Z(t_{j},\boldsymbol{s}_{k})$ は空間データとしてみることができます．
ここで取り上げた例は (データ分析編) でも実際のデータ分析例として取り上げます．

高次元空間データ

時空間データにおいて
$\boldsymbol{Y}(\boldsymbol{s}) = (Z(t_{1},\boldsymbol{s}), \cdots, Z(t_{p},\boldsymbol{s}))'$ とすると(即ち各地点において時間方向の観測をまとめて $p$ 次元のベクトルとみなすと)，時空間データは形式的に

\boldsymbol{Y} = \{\boldsymbol{Y}(\boldsymbol{s}) = (Y_{1}(\boldsymbol{s}),\cdots,Y_{p}(\boldsymbol{s}))' \in \mathbb{R}^{p}: \boldsymbol{s} \in R_{n}\}

と書くことができます(ただし，$Y_{j}(\boldsymbol{s}) = Z(t_{j},\boldsymbol{s})$)．このように変換することで，空間データでも $p \to \infty$ ($n \to \infty$) となることを許せば，空間・時空間データともに $p$ 変量の空間データとして同じ枠組みで扱うことができます．

各地点でのベクトルの次元 $p$ がサンプルサイズ $n$ の増加に従って増加する場合も考えるので，$\boldsymbol{Y}$ は"高次元"の空間データと呼ぶことにします．また観測領域 $R_{n}$ の体積 $\text{vol}(R_{n})$ について $n \to \infty$ として

\text{vol}(R_{n}) = O(\lambda_{n}^{d}), \lambda_{n} \to \infty

の関係が成り立つとしておきます．

空間データに対するブートストラップ法

平均ベクトルの信頼区間の構成

この研究ではデータ $\{\boldsymbol{Y}(\boldsymbol{s}_{i})\}_{i=1}^{n}$ を使ってモデルの平均 $\boldsymbol{\mu} = E[\boldsymbol{Y}(\boldsymbol{s})]$ を推定する問題を考えます. 平均ベクトル $\boldsymbol{\mu} = (\mu_{1},\cdots, \mu_{p})'$ を標本平均

\overline{\boldsymbol{Y}}_{n} = {1 \over n}\sum_{i=1}^{n}\boldsymbol{Y}(\boldsymbol{s}_{i}) =: (\overline{Y}_{1,n},\cdots, \overline{Y}_{p,n})'

で推定するとすると，高次元 ($p = p_{n} \to \infty$) の場合であっても (高次元) 中心極限定理が成り立ち，$n \to \infty$ のとき

\begin{align}
\sqrt{\lambda_{n}^{d}}\overline{\boldsymbol{Y}}_{n} &\approx \boldsymbol{V}_{n} = (V_{1,n},\cdots,V_{p,n})',\\ 
\end{align}

という近似が成り立ちます．ただし，$\boldsymbol{V}_{n}$ は平均ベクトルが $\boldsymbol{0}$ ベクトル，共分散行列が $\Sigma^{\boldsymbol{V}_{n}} = (\Sigma^{\boldsymbol{V}_{n}}_{j,k}) = \lambda_{n}^{d}\text{Var}(\overline{\boldsymbol{Y}}_{n})$ の $p$ 次元の正規分布です¹．ここで，$\text{Var}(\overline{\boldsymbol{Y}}_{n})$ は $\overline{\boldsymbol{Y}}_{n}$ の共分散行列です．以上の結果は標本平均 $\overline{\boldsymbol{Y}}_{n}$ が$\boldsymbol{\mu}=E[\boldsymbol{Y}(\boldsymbol{s})]$ を近似する精度が $\lambda_{n}^{-d/2}(\log p)^{1/2}$，(即ち，$\overline{\boldsymbol{Y}}_{n}-\boldsymbol{\mu}=O_{p}(\lambda_{n}^{-d/2}(\log p)^{1/2})$) であることを意味しています．上記の結果を使うと，さらに $n \to \infty$ で

\sqrt{\lambda_{n}^{d}}\max_{1 \leq j \leq p}\left|{\overline{Y}_{j,n} - \mu_{j}} \over \sqrt{\Sigma^{\boldsymbol{V}_{n}}_{j,j}} \right| \approx \max_{1 \leq j \leq p}\left| V_{j,n} \over \sqrt{\Sigma^{\boldsymbol{V}_{n}}_{j,j}} \right|

という近似が成り立ちます．以下ではこの結果を使って $\boldsymbol{\mu}$ の (同時) 信頼区間を構成することを考えてみたいと思います．まず上記の結果から，$\tau \in (0,1)$ とすると，$\boldsymbol{\mu}$ の $100(1-\tau)$% 信頼区間は

C_{j,n}(1-\tau) = \left[\overline{Y}_{j,n} - \lambda_{n}^{-d/2}\sqrt{\Sigma^{\boldsymbol{V}_{n}}_{j,j}}q_{n}(1-\tau), \overline{Y}_{j,n} + \lambda_{n}^{-d/2}\sqrt{\Sigma^{\boldsymbol{V}_{n}}_{j,j}}q_{n}(1-\tau)\right]

として $C_{n}(1-\tau) = \prod_{j=1}^{p}C_{j,n}(1-\tau)$ の形で与えられます．ここで，$q_{n}(1-\tau)$は

q_{n}(1-\tau) = \inf\left\{x \in \mathbb{R}: P\left(\max_{1 \leq j \leq p}\left| V_{j,n} \over \sqrt{\Sigma^{\boldsymbol{V}_{n}}_{j,j}} \right| \leq x\right) \geq 1-\tau\right\}

で定義される分位点関数です．同時信頼区間 $C_{n}(1-\tau)$において，$\Sigma^{\boldsymbol{V}_{n}}_{j,j}$ や $q_{n}(1-\tau)$ の値は未知なのでこれらをブートストラップ法を用いてデータから推定します．

空間ワイルドブートストラップ

以下ではこの研究で新たに提案した spatially dependent wild bootstrap (SDWB，空間ワイルドブートストラップ) を利用した空間データの平均ベクトルの信頼区間の構成方法について説明します．

(Step1) ブートストラップ標本生成

データ $ \mathcal{Y}_{n} = \{\boldsymbol{Y}(\boldsymbol{s}_{i})\}_{i=1}^{n}$ が与えられたもとで，以下のブートストラップ標本を生成します．

\boldsymbol{Y}^{\ast}(\boldsymbol{s}_{i}) = \overline{\boldsymbol{Y}}_{n} + \left(\boldsymbol{Y}(\boldsymbol{s}_{i} ) - \overline{\boldsymbol{Y}}_{n}\right)W_{n}(\boldsymbol{s}_{i}),\ i = 1,\cdots, n.

ただし，$\{W(\boldsymbol{s}_{i})\}_{i=1}^{n}$ は $E[W_{n}(\boldsymbol{s})] = 0$，$\text{Var}(W_{n}(\boldsymbol{s})) = 1$ かつ

 \text{Cov}(W_{n}(\boldsymbol{s}_{1}), W_{n}(\boldsymbol{s}_{2})) = a\left({\boldsymbol{s}_{1} - \boldsymbol{s}_{2} \over b_{n}}\right)

を満たすような $1$ 変量ガウス過程 $\{W_{n}(\boldsymbol{s}) : \boldsymbol{s} \in \mathbb{R}^{d}\}$ の $\boldsymbol{s}_{1},\cdots, \boldsymbol{s}_{n}$ における離散観測値です．特に $a: \mathbb{R}^{d} \to [0,1]$ は $a(\boldsymbol{0}) = 1$ かつ$[-1,1]^{d}$ にサポートをもつ連続関数で $b_{n}$ は $n \to \infty$ として $b_{n} \to \infty$ を満たす定数列 (データの空間的な従属構造をとらえるバンド幅) です．

(Step2) 共分散行列の推定

このようにブートストラップ標本を生成すると，まず $\Sigma^{\boldsymbol{V}_{n}}$ は

\begin{align}
\hat{\Sigma}^{\boldsymbol{V}_{n}} &= \lambda_{n}^{d}\text{Var}^{\ast}\left(\overline{\boldsymbol{Y}}_{n}^{\ast}\right) \\
&= {1 \over n^{2}\lambda_{n}^{-d}}\sum_{i_{1}, i_{2} = 1}^{n}(\boldsymbol{Y}(\boldsymbol{s}_{i_{1}}) - \overline{\boldsymbol{Y}}_{n})(\boldsymbol{Y}(\boldsymbol{s}_{i_{2}}) - \overline{\boldsymbol{Y}}_{n})'a\left({\boldsymbol{s}_{1} - \boldsymbol{s}_{2} \over b_{n}}\right) \\
\end{align}

で推定できます．ここで $\text{Var}^{\ast}\left(\overline{\boldsymbol{Y}}_{n}^{\ast}\right)$ はデータ $\mathcal{Y}_{n}= \{\boldsymbol{Y}(\boldsymbol{s}_{i})\}_{i=1}^{n}$ が与えられたときのブートストラップ標本平均 $\overline{\boldsymbol{Y}}_{n}^{\ast} = {1 \over n}\sum_{i=1}^{n}\boldsymbol{Y}^{\ast}(\boldsymbol{s}_{i})$ の分散です.

(Step3) 分位点の推定

次に

\overline{\boldsymbol{Y}}^{\ast}_{n} = {1 \over n}\sum_{i=1}^{n}\boldsymbol{Y}^{\ast}(\boldsymbol{s}_{i}) = :(\overline{Y}^{\ast}_{1,n},\cdots,\overline{Y}^{\ast}_{p,n})'

をブートストラップ標本平均とすると，$q_{n}(1-\tau)$ は

\hat{q}_{n}(1-\tau) = \inf\left\{x \in \mathbb{R}: P^{\ast}\left(\sqrt{\lambda_{n}^{d}}\max_{1 \leq j \leq p}\left|{\overline{Y}_{j,n}^{\ast} - \overline{Y}_{j,n}} \over \sqrt{\hat{\Sigma}^{\boldsymbol{V}_{n}}_{j,j}} \right| \leq x\right) \geq 1-\tau\right\}

で推定することができます．ここで，$P^{\ast}$ はデータ $\mathcal{Y}_{n}= \{\boldsymbol{Y}(\boldsymbol{s}_{i})\}_{i=1}^{n}$ が与えられたときのブートストラップ標本の経験分布です²．以上の推定量を $C_{n}(1-\tau)$ に代入すると最終的な $\boldsymbol{\mu}$ の (同時) 信頼区間 $\hat{C}_{n}(1-\tau) = \prod_{j=1}^{p}\hat{C}_{j,n}(1-\tau)$，

\hat{C}_{j,n}(1-\tau) = \left[\overline{Y}_{j,n} - \lambda_{n}^{-d/2}\sqrt{\hat{\Sigma}^{\boldsymbol{V}_{n}}_{j,j}}\hat{q}_{n}(1-\tau), \overline{Y}_{j,n} + \lambda_{n}^{-d/2}\sqrt{\hat{\Sigma}^{\boldsymbol{V}_{n}}_{j,j}}\hat{q}_{n}(1-\tau)\right]

が得られます．このように構成した信頼区間に対して，$n\to\infty$ として

P(\boldsymbol{\mu} \in \hat{C}_{n}(1-\tau)) \to 1-\tau

となることが示せます．

この方法の良い所

(1)高次元の時系列データにも適用可能

[](リモートセンシングの技術進歩により多くの地域で複数の大気汚染物質 (オゾンや $\text{CO}2$，$\text{PM}_{2.5}$ など) の観測が可能になり，観測地点で多くの変数が観測されるようになっています．Nguyen et al.(2017) では各地点において異なる20の気圧層での $\text{CO}2$ 濃度の観測データ ($p=20$) が取り上げられています．このように各地点で数十・場合によっては数百の変数が観測される場合に従来のくうかん)

この研究の設定の良いところは時系列データも特別な場合として扱える点です．空間データの定義で $d = 1$ とすると $\{\boldsymbol{Y}(\boldsymbol{s}) : \boldsymbol{s} \in \mathbb{R}\}$ は "時点" $\boldsymbol{s}_{1},\cdots, \boldsymbol{s}_{n} \in \mathbb{R}$ で $p$ 個の変数を観測する時系列データとしてみることができるので，各時点において観測される変数の数 $p$ が時系列の観測期間 $n$ よりも大きいような 高次元時系列データにも対応できるところです．高次元時系列データの例としては金融データがあります．金融データの場合，各時点で多数の金融資産価格を同時に観測する場合に対応します．

(2)時空間データに適用する場合，時間方向の制約がほとんどなく，様々なモデルに対応可能

本研究では時空間データに対して，空間方向には定常な空間モデルで漸近的に近似できることを仮定しますが，時間方向については非定常の離散時間 or 連続時間モデルであっても空間ワイルドブートストラップを用いたデータ分析が可能で，理論的にもその妥当性を示すことができます．

(3)実装が容易

空間データに対する既存のブートストラップ法と比べて，今回提案した空間ワイルドブートストラップ法はチューニングパラメータがガウス過程のバンド幅 $b_n$ のみで，より実装が簡単になっています³.
以下でより具体的なアルゴリズムについてもまとめておきます：
(Step1) 標本平均 $\overline{\boldsymbol{Y}}_{n}$ と $\hat{\Sigma}^{\boldsymbol{V}_{n}}$ を計算．
(Step2) ブートストラップ標本の数 $B$ を設定し，共分散行列 $A = (a(\boldsymbol{s}_{i}-\boldsymbol{s}_{k})/b_{n})_{1 \leq i,k \leq n}$ をもつ平均 $\boldsymbol{0}$ の $n$ 次元正規分布に従う変数 $\boldsymbol{W}^{(1)}$$,\dots,$ $\boldsymbol{W}^{(B)}$ を生成：

\boldsymbol{W}^{(k)} = (W_{1}^{(k)},\dots,W_{n}^{(k)})',\ k=1,\dots, B.

(Step3) ブートストラップ標本 $\mathcal{Y}^{(1)}_{n}, \dots, \mathcal{Y}^{(B)}_{n}$ を生成：

\begin{align}
\mathcal{Y}^{(k)}_{n} &= \{\boldsymbol{Y}^{(k,\ast)}(\boldsymbol{s}_{1}),\dots, \boldsymbol{Y}^{(k,\ast)}(\boldsymbol{s}_{n})\},\\
\boldsymbol{Y}^{(k,\ast)}(\boldsymbol{s}_{i}) &= \overline{\boldsymbol{Y}}_{n} + \left(\boldsymbol{Y}(\boldsymbol{s}_{i} ) - \overline{\boldsymbol{Y}}_{n}\right)W^{(k)}_{i},\ k=1,\dots,B.
\end{align}

(Step4) ブートストラップ標本 $\mathcal{Y}_{n}^{(k)}$ ごとに

T_{n}^{(k)} = \sqrt{\lambda_{n}^{d}}\max_{1 \leq j \leq p}\left|{\overline{Y}_{j,n}^{(k,\ast)} - \overline{Y}_{j,n} \over \sqrt{\hat{\Sigma}_{j,j}^{\boldsymbol{V}_{n}}}}\right|,\ k=1,\dots, B

を計算して $T_{n}^{(1)},\dots, T_{n}^{(B)}$ の経験分布から分位点 $q_{n}(1-\tau)$ の推定値 $\hat{q}_{n}(1-\tau)$ を計算．

(4)観測地点が格子上でも不等間隔でも実装方法が変わらない

ブートストラップ標本の生成にあたり，(i)標本平均の計算と(ii)ガウス過程を使って，各地点でデータをランダムに発生させる，という手順をとるので観測地点がどのように配置されていても実装方法は変わりません．

より詳しい設定や理論的背景については論文をご参照ください．

まとめ

この記事では「空間データのためのブートストラップ法」について解説しました．具体的には空間ワイルドブートストラップ法とその空間データ・時空間データへの適用方法について主に理論的側面から紹介しました．実際のデータ分析の例については別の記事（データ分析編）で紹介したいと思います．
株式会社Nospareには空間統計・時空間統計に限らず，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospare までお問い合わせください．

参考文献
[1] Kurisu, D., Kato, K. and Shao, X. (2021) Gaussian approximation and spatially dependent wild bootstrap for high-dimensional spatial data. Journal of the American Statistical Association 119, 1820-1832. arXiv:2103.10720.

(補足) 正確には $\boldsymbol{V}_{n}$ の分布は観測地点に依存するので，観測地点 $\boldsymbol{S}_{n} = \{\boldsymbol{s}_{1},\cdots, \boldsymbol{s}_{n}\}$ が与えられたもとで $\boldsymbol{V}_{n}$ は (観測地点に依存するランダムな分散共分散行列を持つ) 正規分布になります． ↩
(補足) ここでもより正確には観測地点の集合 $\boldsymbol{S}_{n}$ とデータ $\mathcal{Y}_{n}$ が与えられたときの条件付分布を考える必要があります．実用上は特に気にする必要はありませんが，数学的に正確に記述しようとすると複雑になってしまいます．．． ↩
(補足) 時系列データや空間データのように従属性をもつデータに対するブートストラップ法としてはこの研究で提案した方法 (wild bootstrap あるいは multiplier bootstrap と呼ばれます) 以外にも block bootstrap と呼ばれる方法があります．後者の方法はデータをある範囲ごとにまとめて "block" を構成してそれらの block をリサンプリングする方法で，リサンプルされた block 内ではデータの従属性が保たれるのでもともとのデータの従属性を損なわずに推定量の分位点が推定することができるようになります．Wild bootstrap では block bootstrap と比べ，データを block に分割する必要ないのでチューニングも少なく(バンド幅 $b_n$ を決めるだけ)，より実装が簡単になります． ↩