バリオグラム

概要

アメダスデータを基にした地図の変形プログラムでは,観測地点間の観測値の相関の大小をもとに,地点間の距離を変化させた.この際に注意が必要なのは,距離が近い観測点間は当然観測値の相関が高いため,2点間の「距離の割には」相関が高いもしくは低いという規準が必要となる.そのためには,2点間の距離に応じたもっともらしい相関の値を推定する必要があるだろう.そこで有効なのが,空間統計学という分野で用いられてきたバリオグラム(Variogram)である.ここでは,空間的な相関を推定するためのバリオグラムの理論について説明する.

バリオグラムとは

平面上の各点\(\boldsymbol{x}\in\mathbb{R}^2\)に,確率的な値\(Y(\boldsymbol{x})\in\mathbb{R}\)が割り当てられているとしよう(こういったものを確率場という).気象データの例でいうと,\(\boldsymbol{x}\)は必ずしも観測地点とは限らない日本地図上の各点であり,\(Y\)はそこでの降雨量などの気象値である.つまり,観測地点以外での気象値も含めて統計モデルを構成する.このとき,バリオグラム\(\gamma:\mathbb{R}^2\times \mathbb{R}^2\rightarrow \mathbb{R}_{\geq 0}\)は二点間の\(Y\)の値の差の分散として定義される. \[ \gamma(\boldsymbol{x},\boldsymbol{x}')={\rm Var}(Y(\boldsymbol{x})-Y(\boldsymbol{x}')). \] また,分散を分解すると以下のように書き換えることもできる. \[ \gamma(\boldsymbol{x},\boldsymbol{x}')={\rm Var}(Y(\boldsymbol{x}))+{\rm Var}(Y(\boldsymbol{x}'))-2{\rm Cov}(Y(\boldsymbol{x}),Y(\boldsymbol{x}')). \]

特に\(Y(\cdot)\)が定常的,つまり任意の\(\boldsymbol{x},\boldsymbol{h}\in\mathbb{R}^2\)に対して\(Y(\boldsymbol{x})\)と\(Y(\boldsymbol{x+h})\)の確率分布が一致するとき,\(\boldsymbol{h}=\boldsymbol{x}-\boldsymbol{x}'\)とすると,\(\gamma(\boldsymbol{x},\boldsymbol{x}')=\gamma(\boldsymbol{x}-\boldsymbol{x},\boldsymbol{x}'-\boldsymbol{x})=\gamma(0,\boldsymbol{h})\)となり,これは\(\boldsymbol{h}\)のみに依存する.やや表記の乱用であるが,これを改めて\(\gamma(\boldsymbol{h})\)と書くことにする.また,各点での分散\({\rm Var}(Y(\boldsymbol{x}))\)は\(\boldsymbol{x}\)に依存しなくなるため,これを単に\(V\)と表すと以下のようにバリオグラムと相関係数の関係が求まる.

バリオグラムと相関係数との関係(定常性の過程のもと) \[ \gamma(\boldsymbol{x},\boldsymbol{x}')=2V-2{\rm Cov}(Y(\boldsymbol{x}),Y(\boldsymbol{x}'))=2V\{1-{\rm Cor}(Y(\boldsymbol{x}),Y(\boldsymbol{x}'))\}. \] \[ {\rm Cor}(Y(\boldsymbol{x}),Y(\boldsymbol{x}'))=1-\frac{\gamma(\boldsymbol{x},\boldsymbol{x}')}{2V}. \tag{1} \]

さらに,\(Y(\cdot)\)が等方的,つまり任意の\(\boldsymbol{x},\boldsymbol{h}\in\mathbb{R}^2\)に対して\(y(\boldsymbol{x})-y(\boldsymbol{x+h})\)の確率分布が,\(\boldsymbol{h}\)の方向にはよらずノルムにのみ依存すると仮定する.このとき,\(\gamma(\boldsymbol{x},\boldsymbol{x}+\boldsymbol{h})=\gamma(0,\boldsymbol{h})\)を\(\boldsymbol{h}\)のノルム\(r=\|\boldsymbol{h}\|\)を用いて\(\gamma(r)\)と書くことにする(表記の乱用2回目).気象データの解析ではこのように定常性と等方性を仮定したもとで,\(\gamma(r)\)と相関係数の関数を推定した.

バリオグラムの理論

バリオグラムはその定義より非負の値をとり,さらに\(\gamma(0)=0\)をみたす.しかし,このようなすべての関数が何らかの確率場\(Y(\cdot)\)のバリオドグラムとなり得るわけではない.実は,関数がバリオクラムとなるための必要十分条件として,以下のことが知られている.

バリオグラムの条件付き負定値性 \(\gamma(\cdot,\cdot)\)がバリオグラムとなるための必要十分条件は,条件付き負定値対称関数であることである.つまり,\(\gamma(\boldsymbol{x},\boldsymbol{x}')=\gamma(\boldsymbol{x}',\boldsymbol{x})\)であり,かつ有限個の実数の組でその和が0になるような任意の\(z_1,\dots,z_N\)および,任意の\(x_1,\dots,x_N\in \mathbb{R}^2\)に対して以下が成立する: \[ \sum_{i=1}^N \sum_{j=1}^N z_i z_j\gamma(\boldsymbol{x}_i,\boldsymbol{x}_j)\leq 0. \]

条件付き負定値行列は,その名の通り\(\sum_{i=1}^N z_i=0\)という制約を抜かすと負定値行列(正定値行列のマイナス)の定義と一致する.つまり,負定値行列であれば条件付き負定値行列である(ふつうは条件を付けると狭くなるが,定義内の仮定に条件を付けているので逆になることに注意).条件付き負定値行列は正定値行列と同様に,カーネル法などの再生核ヒルベルト空間での統計学で重要な役割を果たすことが知られている.実際,バリオグラムの理論をはじめとする空間統計学のkrigingと呼ばれる手法は,再生核ヒルベルト空間理論の統計学への初期の応用分野のひとつである.

指数バリオグラム

定常性と等方性を仮定すると,バリオグラムは\(\gamma(r)\)のように書けるのであった.一方,指数関数\(\exp(-r)\)は正定値関数であるため,正の数\(c\)を用いた \[ \gamma(r)=1-\exp(-cr) \] は負定値関数で,バリオグラムとなる.一方,\(r=0\)でのみ\(0\)で,\(r>0\)で\(1\)をとる関数\(\boldsymbol{1}_{(0,\infty)}(r)\)も負定値関数となることがすぐわかり,これもバリオグラムになる.バリオグラムの条件付き負定値性を用いた別表現からもわかるように,バリオグラムを正のスカラー倍したり,また二つのバリオグラムを足し合わせた結果もまたバリオグラムとなる.よって,正の値\(s,n\)を用いた以下の関数も負定値関数であり,バリオグラムとなる.

指数バリオグラム \[ \gamma(r)=(s-n)\{1-\exp(-cr)\}+n\boldsymbol{1}_{(0,\infty)}(r) \tag{2} \]

ここで,\(n\)は\(\lim\downarrow 0\)の極限でのバリオグラムの値で,ナゲット(nugget)と呼ばれる.一方,\(s\)は\(\lim\rightarrow \infty\)の極限でのバリオグラムの値で,シル(sill)と呼ばれる.このようにバリオグラムは原点でジャンプを持ってもよいが,条件付き負定値関数の議論からジャンプが存在するのは原点のみであり,原点以外では連続であることも証明できる.なお,ここでは指数バリオグラムのみ説明したが,球面的バリオグラムやガウシアンバリオグラムなど多種のバリオグラムが応用上知られいている.

バリオグラムの地図変形への応用

アメダスデータの地図変形プログラムでは,以下のような回帰モデルで相関と距離との関係式を推定した. \[ \begin{aligned} r(i,j) &= a + b\exp(-c \cdot d(i,j)) + \varepsilon_{(i,j)} \\ \varepsilon_{(i,j)} &\sim {\rm N}(0, \sigma^2) \end{aligned} \] この回帰式は,定常性,等方性の仮定のもとでの指数バリオグラムより導かれる相関の式と一致する.実際,式(1),(2)より \[ \begin{aligned} {\rm Cor}(Y(\boldsymbol{x}),Y(\boldsymbol{x}'))&=1-\frac{\gamma(r)}{2V} \\ &= 1 -\frac{s-n}{2V}\{1-\exp(-cr)\}-\frac{n}{2V}\boldsymbol{1}_{(0,\infty)}(r) \\ &=   \begin{cases} 1 & (r=0),\\ \frac{1}{2V}\{(1-s)+(s-n)\exp(-cr)\} & (r>0). \end{cases} \end{aligned} \] となり,さらに \[ {\rm Cor}(Y(\boldsymbol{x}),Y(\boldsymbol{x}')) \mapsto r(i,j)\] \[ r \mapsto d(i,j) \] \[ \frac{1-s}{2V} \mapsto a \] \[ \frac{s-n}{2V} \mapsto b \] のように置きなおすと,\(r(i,j) = a + b\exp(-c \cdot d(i,j))\)の回帰式の形となることがわかる.
なお,地図変形プログラム用のデータ解析では二乗誤差最小化で推定を行ったため,モデルのノイズは正規分布としているが,より正確な統計解析が必要なときは非正規分布のノイズについても候補として考えたほうが良いことを付記しておく.