« 読了：Rabiee, F. (2004) グループ・インタビューの逐語録をどうやって分析するか | メイン | 読了：Bowman & McMahan (2007) ヴァーチャル・リアリティってのは視覚的没入性が高ければ高いほどよいというものでもないだろう »
2017年9月14日 (木)
Zammit-Mangion, A., Cressie, N. (2017) Fixed Rank Kriging: The R Package.
Rの空間統計パッケージFRKのvignette。適当にめくるだけのつもりが、面白くてついつい読み耽ってしまった。ま、最初の3頁だけだけど。
FRKとはFixed Rank Kriging (固定ランク・クリギング)の略である。
いわく。
FRKと似たパッケージとして以下がある。[このくだり、まったく意味がわからない箇所もあるが、後日のために逐語訳する。精度行列というのは共分散行列の逆行列のことね]
- LatticeKrig: Wendland 基底関数(コンパクトなサポートを持つ)を使って、空間的相関を持つ過程を分解する。マルコフ想定によって精度行列(後述する$\mathbf{K}^{-1}$)を構築し、これらの基底関数の係数のあいだの依存性を記述する。我々が言うところの微細スケールプロセス変動[BAU内変動のこと?]を提供するのではなく、過程のスケールの微細さは、利用されている基底関数の微細さが限度となる。しかし、計算上のモチベーションにより$\mathbf{K}^{-1}$にスパース性が付与されているので、このスケールは比較的に微細なものとなりうる。基盤となるモデルは、ガウシアン・マルコフ確率場(GMRF)仮定を用いて構築されたスパースな精度行列を用いており、その結果、計算は効率的であり、また多数(10000以上)の基底関数が使えるようになっている。
- INLA: モデル・フィッティングと予測のための汎用パッケージ。空間モデリング・時空間モデリングに使用する場合は次のように仮定するのが一般的である。まず、基底関数として三角形の「テント」関数を仮定する。係数は正規分布しその精度行列はスパースだと仮定する。その結果としてガウシアン過程の共分散関数はMaternクラスの空間共分散関数の近似となる。このように、INLAのアプローチはLatticeKrigと多くの点で共通している。INLAの主要な利点は、いったん空間モデルないし時空間モデルをつくってしまえば、パッケージが提供しているすべての近似推論機構と尤度モデルが使えるようになる、という点である。
- spBayes: Kang & Cressie (2011)が開発したベイジアンFRKでは、空間基底関数は固定され$\mathbf{K}$に事前分布が付与される。Banerjee et al (2008)の予測過程アプローチもベイジアンFRKの一種だとみなすことができる。このアプローチでは、基底関数は空間ランダム効果の共分散関数から構築され、従ってパラメータに依存する。この予測過程を実装したのがspBayesパッケージである。多変量の空間・時空間過程を扱うことができる。多様な尤度モデルについて、MCMCによるベイズ推論を行う。結果として得られる基底関数はパラメトリックな共分散モデルに基づいて構築されているので、パラメータについての事前分布に基づき、MCMCの反復ごとに、新しい基底関数が生成される。そのため計算が遅くなることがあり、予測過程におけるノットの数を小さくしておく必要がある。
モデルの概要。
FRKパッケージはふつうのバリオグラム・モデルではなく、空間ランダム効果モデル(SRE)で問題を定式化する。
地点$\mathbf{s}$における量$Y(\mathbf{s})$に関心があるとしよう。古典的な空間モデルなら、共変量ベクトルを$\mathbf{t(s)}$、回帰係数ベクトルを$\mathbf{\alpha}$、空間的相関があるランダム効果を$\upsilon(\mathbf{s})$、空間的相関がないランダム効果を$\xi(\mathbf{s})$ととして
$Y(s) = \mathbf{t(s)}^\mathrm{T} \mathbf{\alpha} + \upsilon(\mathbf{s}) + \xi(\mathbf{s})$
ただし$E(\upsilon(\cdot)) = E(\xi(\cdot)) = 0$、とするところである。
ここで、
$\upsilon(\mathbf{s}) = \sum_{l=1}^r \phi_l(\mathbf{s}) \eta_l$
とする。$\mathbf{\eta} \equiv (\eta_1, \ldots, \eta_r)^{\mathrm{T}}$はランダムベクトル、$\mathbf{\phi} \equiv (\phi_1(\ldots), \ldots, \phi_r(\ldots))$はあらかじめ定められた空間基底関数である。
[理解するまでに手間取ったのだが、ここがこの話のミソなのである。この$\mathbf{\phi}$こそがバリオグラム・モデルの変わり果てたお姿なのだと思う。話の先取りになるけど、FRKパッケージにはデータとバリオグラム関数の形状を指定すると$\mathbf{\phi}$の行列を自動的に作ってくれる関数がある。それを使わずに自力で作ってもよい]
さて、空間を$N$個の重ならない小さな空間に完全に分ける[メッシュみたいなものであろう]。これを基本地域単位 BAU と呼ぶ。BAUの数$N$は基底の数$r$よりずっと大きいものとする。
それぞれのBAUにおいて$Y(\mathbf{s})$を平均する。つまり、$i$番目のBAU $A_i$において
$Y_i \equiv \frac{1}{|A_i|} \int_{A_i} Y(\mathbf{s}) d(\mathbf{s})$
このレベルでも、
$Y_i = \mathbf{t}_i^\mathrm{T} \mathbf{\alpha} + \upsilon_i + \xi_i$
と分解できる。
以下、各項について順にみていくと...
BAUレベルの共変量$\mathbf{t}_i$はこうなる。
$\mathbf{t}_i \equiv \frac{1}{|A_i|} \int_{A_i} \mathbf{t}(\mathbf{s}) d \mathbf{s}$
BAUレベルの相関ありランダム効果$\upsilon_i$は
$\upsilon_i \equiv \frac{1}{|A_i|} \int_{A_i} \upsilon(\mathbf{s}) d \mathbf{s}$
ここに$\upsilon(\mathbf{s}) = \mathbf{\phi}(\mathbf{s}) \mathbf{\eta}$を代入すると... [面倒くさいから書かないけど] 係数ベクトルと$\mathbf{\eta}$の積になる。係数の部分を取り出して$N \times r$行列$\mathbf{S}$とすれば、$\mathbf{\upsilon} = \mathbf{S\eta}$となる。なおこのパッケージでは、実際にはそれぞれの係数について積分するのではなく、BAUは十分に小さいとみて、BAUの重心$\mathbf{s}_i$を使い$\mathbf{\phi}(\mathbf{s}_i)$と近似する。
$\mathbf{\eta}$について。このパッケージでは、$\mathbf{\eta}$を平均0, 共分散行列$\mathbf{K}$の正規ベクトルとする。$\mathbf{K}$は自由推定してもいいし、なんらかの構造を与えても良い。前者をFRK-Vと呼び(Vはバニラの略)、後者をFRK-Mと呼ぶ(Mはモデルの略)。
相関なしのランダム効果について。$\xi(\mathbf{s})$についてはもう関心を持たない。BAUのレベルでの平均$\xi_i$について、平均0, 分散$\sigma^2_\xi \nu_{\xi i}$の正規分布に独立に従うとする。ここで$\sigma^2_\xi$は全BAUを通したパラメータ、$\nu_{\xi i}$は既知の定数でBAUの異質性を表す。
ここまでをまとめよう。BAUレベルの式
$Y_i = \mathbf{t}_i^\mathrm{T} \mathbf{\alpha} + \upsilon_i + \xi_i$
を書き直して
$\mathbf{Y} = \mathbf{T} \mathbf{\alpha} + \mathbf{S} \mathbf{\eta} + \mathbf{\xi}$
$\mathbf{T}$は行数$N$の共変量行列、$\mathbf{\alpha}$は係数ベクトル。$\mathbf{S}$は$N \times r$の行列、$\mathbf{\eta}$は長さ$r$のランダム効果ベクトルでその共分散行列が$\mathbf{K}$。$\mathbf{\xi}$は長さ$N$のランダム効果ベクトルでその共分散行列は$\sigma^2_\xi \mathbf{V}_\xi$、ただし$\mathbf{V}_\xi$は既知。
さて、ここまで考えてきた$Y(\cdot)$というのは潜在過程であって、観察自体は$m$個のフットプリントについてのみ可能であるとしよう。ここでフットプリントとは、ひとつ以上のBAUからなる地域のことで、重複していてもかまわない。$m$は$r$よりずっと大きいものとする。なお、$m$はBAUの数$N$より大きくても小さくてもよい。
フットプリント$B_j$における観察値は次の3つの和とする[面倒くさいので式は省略]：
- $B_j$にそこに含まれるBAUの$Y_i$の平均。
- $B_j$に含まれるBAUのそれぞれが持つ体系的誤差を表す$\delta_i$の平均。平均0, 分散$\sigma^2_\delta \nu_{\delta i}$の正規分布に独立に従うものとする。ここで$\nu_{\delta i}$は既知。
- $B_j$の測定誤差$\epsilon_j$。平均0の正規分布に独立に従うものとする。[この分散はパラメータにしない模様]
結局、推定するパラメータは、バニラFRKでは$\mathbf{\alpha}, \sigma^2_\xi, \sigma^2_\delta$, そして$\mathbf{K}$ の4つ。モデルFRKでは、$\mathbf{K}$の代わりに$\mathbf{K}$について組んだモデルのパラメータがはいることになる。
[以下、力尽きたので、ほとんど読んでいない]
パラメータはEMアルゴリズムで推定する...
このモデルにより、任意の予測領域についての予測が可能で...
コード付きの事例が2つ...
時空間でクリギングする事例が2つ...
空間異方性を持たせるには...
既定関数とBAUを手動で与えるには...
今後の課題:
- 現状ではFRKはanalytic formによるlocal basis functionsを用いているのだが、関数形が未知の基底関数、たとえば経験直交関数などにも対応したい...[ごめんなさい、なにいってんだかさっぱりわかんない]
- 現状ではBAUより下の分散は無視してるんだけど、これを拡張して...
- 三次元超平面とかに拡張して...
- 現状では地点数が数十万になると速度が落ちるので...
- BAUの大きさが違う場合に拡張して...
云々、云々。
...というわけで、正直なところ頑張って読んだのは最初の3頁だけなんだけど、固定ランク・クリギングという発想がよくわかった(ような気がする)。$m$個の観察値の後ろに$N$個のメッシュの値を考え、そのまた後ろにもっと少ない$r$個の隠れた値を考えるわけね。時系列モデルで言うと、観察変数のうしろに状態変数を考え、そのまた後ろにたまにしか動かないような状態変数を考えているわけだ。
だから大データでもクリギングが楽だ、ってわけね。なるほどねー、面白いなあ。(←すっかりわかったような気分)
論文：データ解析(2015-) - 読了：Zammit-Mangion & Cressie (2017) 大きな空間データを固定ランク・クリギングするRパッケージ FRK