特集「時空間統計解析：新たなる分野横断的展開」［総合報告］

統計数理（2012）
第 60 巻 第 1 号 57–71
c 2012 統計数理研究所


時空間大規模データに対する統計的解析法†

矢島 美寛 ・平野 敏弘†（受付 2011 年 7 月 1 日；改訂 9 月 1 日；採択 9 月 8 日）

要旨

本論文では，時空間大規模データを効率的に解析するために提案された統計的解析法の最近の発展についてサーベイする．まず最良不偏線形予測量を近似し，迅速に予測量を計算するために用いられる Covariance Tapering について解説する．次に時空間領域および周波数領域において最尤推定法を代替し，高速計算可能な推定法を紹介する．次にパラメータ数を縮約し少数の因子によりデータを説明する潜在過程モデルを紹介する．最後に今後解決すべき課題を列挙する．
キーワード： 時空間統計解析，大規模時空間データ，Covariance Tapering，潜在過程
モデル．
1.

序

近年時空間データを収集・整理・解析するためのインフラストラクチャーはリモートセンシ
ング，全地球測位システム（Global Positioning System, GPS）
，地理情報システム（Geographical
Information System, GIS）などの科学技術の発達，また R や Matlab など統計解析ソフトウエ
アの普及により整備され充実してきた．
しかしながら時空間データの総数は千，万のオーダーになることは頻繁であり，さらには数十万から百万のオーダーに上ることも稀ではない．したがって大容量の記憶能力および高速の処
理能力を保持する現在の計算機をもってしても，このようなサンプル数のデータを迅速に解析することが困難になる場合が生じる．たとえばパラメトリック・モデルの推定や最良線形予測量の構成には，与えられた観測値間の共分散行列の逆行列を計算する必要がある．通常のアルゴリズムでは逆行列の計算に要するオペレーション回数は，サンプル数の 3 乗に比例する．したがってサンプル数が 1,000 個でもオペレーションの回数は 10 億のオーダーになってしまう．
本論文では，21 世紀に入りこの 10 年ほどの間に発展してきた時空間大規模データの効率的な統計的解析法について報告する．
ここで以降の章の構成を述べておく．2 章では準備として時空間データの数学的表現（確率場と言う），定常確率場およびそれに対する代表的なモデルを説明する． 3 章では定常確率場の予測やパラメータ推定に用いられる Covariance Tapering について解説する．4 章では不規則間隔で観測される時空間データに対して，最尤推定法に代わる簡便で効率的な推定法のいくつかを紹介する．5 章ではパラメータ数を極力減らして，データ解析を迅速にする目的で提案された潜在過程モデル（Latent process model）を紹介する．最後に 6 章では今後の課題を列挙する．
なお関連したトピックに関するコンパクトかつ明解なサーベイ論文として Su et al.（2011）が†東京大学大学院 経済学研究科：〒113–0033 東京都文京区本郷 7–3–1

統計数理

58

第 60 巻

第 1 号 2012

ある．
2.

時空間データの数学的表現・確率場

本章ではまず時空間データをどのように数学的に表現するかについて説明する．次にこの表現に基づいて定常確率場及びそれに対する代表的なモデルを紹介する．
実数の全体を R = (−∞, ∞) とし，その d(= 1, 2, . . . ) 次元ユークリッド空間は Rd と表す．また整数の全体を Z = {0, ±1, ±2, . . . } とし，その d 次元直積集合 Z × · · · × Z は Z d と表す．両者を統一的に表す場合には K d とする．
次に観測地点（site）
を s(∈ K d ) とし，そこで観測されるデータは確率変数 Y (s) と表す．Y (s)がスカラーのときは一変量データ，ベクトルのときは多変量データである．本論文では一変量データを扱う．例えば d = 2 のときには，s は 2 次元ベクトルで第 1，2 座標は各々緯度，経度を示し，Y (s) はその地点における地価とする．また d = 3 のときには，s は 3 次元ベクトルで，第 1，2，3 座標は各々緯度，経度，高さを示し，Y (s) はその地点における気温などとする．sの動く領域を D(⊂ K d ) としたとき，データの全体を {Y (s) : s ∈ D} と書き，確率場（random
と言う．
ﬁeld）
さらに観測時点も考慮する場合には，次元 d を一つ大きくし s の最後の座標が時点を表すとすればよい．先ほどの気温の例では d = 4 として，s の第 4 座標を時点に取ればよい．時点を強調したいときには第 4 座標のみ分離して t(∈ K) と記し，データを Y (s, t) と表す．
確率場 {Y (s) : s ∈ K d } が次の 2 条件をみたすとき，弱定常確率場（weakly stationary random
ﬁeld）
と言う．
（i）期待値は s に依存せず一定の値 E(Y (s)) = µ を取る．以下では簡単のため断りのない限
り µ = 0 とする．
（ii）共分散はベクトル差 t − s のみに依存し，C(t − s) = E(Y (t)Y (s)), t, s ∈ K d となる．
{C(h), h ∈ K d } を自己共分散関数（autocovariance function）と言う．以下では弱定常確率場
を単に定常確率場と呼ぶ．特に d = 1 の場合は時系列解析における定常過程になる． Y (s) および C(h) はスペクトル表現が可能で，それぞれY (s) =
exp(is λ)dM (λ)
Td
C(h) =(2.1)
exp(ih λ)dF (λ)
Tdと表すことが出来る．ここで i は虚数単位 i2 = −1 であり， s = (s1 , . . . , sd ) ， h = (h1 , . . . , hd ) ，
λ = (λ1 , . . . , λd ) とし， は転置ベクトルを意味する．K = Z のときは T = (−π, π]，K = R のときは T = (−∞, ∞) である．M = {M (λ), λ ∈ T d } は d 次元複素直交増分過程， F (λ) は d 次
元非負測度でスペクトル分布関数（spectral distribution function）と各々呼ばれ，任意のボレル
集合 ∆, ∆1 , ∆2 (∆1 ∩ ∆2 = φ) ⊂ T d に対して，E|M (∆)|2 = F (∆), E(M (∆1 )M (∆2 )) = 0 が成立
する．スペクトル分布関数が絶対連続な場合には，その密度関数 f (λ) をスペクトル密度関数（spectral density function）と言う（Guyon, 1995; Stein, 1999; Yaglom, 1987）．このとき（2.1）は(2.2)
C(h) =
exp(ih λ)f (λ)dλTdになる．
定常確率場に対して統計モデルを導入する方法は大別して二つに分かれる．ひとつは Y (s)に対して直接導入する方法であり，他方は自己共分散関数 C(h) あるいはスペクトル密度関数

時空間大規模データに対する統計的解析法

59

f (λ) に導入する方法である．前者にも有用なモデルがあるが，d ≥ 2 の場合になると，d = 1 の場合すなわち定常時系列のように過去・現在・未来と言う自然な順序がなく，一般にはその意
味づけあるいは数学的正当性の証明が難しい．一方自己共分散関数は非負定値関数であればよく，その必要十分条件はスペクトル表現（2.1）
を持つことである．またその形状から地点間の観測値の相関関係も把握しやすい．したがってここでは後者のアプローチを採用する．
現時点までに提案されているモデルの中で代表的なものとして等方型モデル（Isotropic Model）と分離型あるいは乗法型モデル（Separable Model）がある．
d2
等方型モデルは，自己共分散関数 C(h) が距離 ||h|| =i=1 hi のみに依存し，方向には無関係なモデルである．すなわち C0 (x)(x ∈ R) を 1 変数正定値関数としたとき，C(h) は x = ||h||を代入してC(h) = C0 (||h||)によって表現される．
C0 (x) に対するパラメトリック・モデルのなかで，応用上頻繁に用いられるモデルは Matérn族（Matérn, 1947, 1960, 1986; Cressie, 1993; Stein, 1999）である．Matérn 族に関する歴史については Guttorp and Gneiting（2006）を参照されたい．Matérn 族の一般形は
C0 (x) =π 1/2 φ2ν−1 Γ(ν + 1/2)α2ν

(α|x|)ν Kν (α|x|)によって表現される．ここで α, ν, φ は正の定数，また Kν は変形された Bessel 関数（modiﬁed
Bessel function）と呼ばれる特殊関数である．φ が共分散の大きさ，α が形状，ν が滑らかさお
よび x → ∞ のときの 0 への収束速度を規定するパラメータである．例えば ν = 1/2 のときは
K1/2 (x) = π/(2x) exp(−x)， ν = 3/2 のときは K3/2 (x) = π/(2x)(1 + x−1 ) exp(−x) になる．
したがって自己共分散関数は ν = 1/2 のときは C0 (x) = πφα−1 exp(−α|x|)，ν = 3/2 のときは
C0 (x) = 12 πφα−3 exp(−α|x|)(1 + α|x|) になる．この自己共分散関数に対応するスペクトル密度関数は ω = ||λ|| のみに依存しf (ω) =φπ (d−1)/2 (α2 + ω 2 )ν+d/2

によって表現される．
次に分離型モデルでは，C(h) が 2 個以上の正定値関数の積として表現される．例えば h を
2 つのベクトル h̃1 = (h1 , . . . , hm ) , h̃2 = (hm+1 , . . . , hd ) に分割し，C(h) は 2 つの正定値関数
C1 (h̃1 ), C2 (h̃2 ) の積 C(h) = C1 (h̃1 )C2 (h̃2 ) によって定義される．
この統計モデルを周波数領域から見ると，λ = (λ1 , . . . , λd ) を 2 つのベクトル λ̃1 = (λ1 , . . . , λm ) ,λ̃2 = (λm+1 , . . . , λd ) に分割し，スペクトル密度関数 f (λ) は 2 つの正の値を取る関数 f1 (λ̃1 ), f2 (λ̃2 )の積 f (λ) = f1 (λ̃1 )f2 (λ̃2 ) によって表現される．
ただし等方型あるいは分離型モデルが適さない現実のデータは多々ある．現在非等方型・非分離型のモデルの開発も精力的に行われていることを断っておく．関心のある読者は矢島（2011）
およびそれに所収の参考文献を参照されたい．
3.

Covariance Tapering の予測への応用

Covariance Tapering とは，真の自己共分散関数をサポートが有限な自己共分散関数で代替
し，さまざまな統計的推測を行う方法である．時系列解析とのアナロジーで言えば，真のモデ
ルをある一定の時間間隔以上離れた観測値が互いに無相関になる移動平均モデルで近似して統計的解析を行うことに対応する．本章では Covariance Tapering の応用例として既知の観測値

60

統計数理

第 60 巻

第 1 号 2012から未知の値を予測する問題を考える．なお推定への応用については紙幅の関係で割愛するが，
Du et al.（2009）, Kaufman et al.（2008）, Wang and Loh（2011）などの文献を参照されたい．
まず最良線形不偏予測量（Best Linear Unbiased Predictor，BLUP）を定式化する．次に特定化の誤りによりあるいは意図的に，真の自己共分散関数とは異なった自己共分散関数を用いて構成された BLUP の最小 2 乗誤差と真の BLUP のそれが漸近的に等しくなるための十分条件を説明する．次に Covariance Tapering について説明し，この方法に基づく BLUP が十分条件を満たす場合を明らかにする．
いま回帰モデルY (s) = m(s) β + (s)を考える．ここで m(s)，β は p 次元ベクトルで各々説明変数，回帰係数ベクトルとし，{(s)}は期待値 0 の定常確率場とする．このとき Y = (Y (s1 ), Y (s2 ), . . . , Y (sn )) が観測され，別の地
点 s∗ の値 Y (s∗ ) を予測するとしよう．観測値に対して線形 λ Y (λ ∈ Rn ) でありかつ不偏性E(λ Y ) = λ M β = m(s∗ ) β をみたす予測量の中で平均 2 乗誤差 E[Y (x∗ ) − λ Y ]2 を最小にする
予測量を BLUP と言う．ここで M は n × p 行列 M = (m(s1 ), . . . , m(sn )) とする．
いま n × n 行列 C ，n 次元ベクトル c∗ をC = (cij ), cij = Cov((si ), (sj ))c∗ = (c∗1 , . . . , c∗n ) , c∗i = Cov((s∗ ), (si ))によって定義する． C は正則行列，M はフル・ランク行列すなわち rank(M ) = p のとき，BLUP
およびその平均 2 乗誤差は(3.1)

Ŷ (s∗ ) = (c∗ ) C −1 (Y − M β̂) + m(s∗ ) β̂

(3.2)

E[Y (s∗ ) − Ŷ (s∗ )]2 = c0 − (c∗ ) C −1 c∗ + γ  M  C −1 M γ

によって与えられる．ここで β̂ は回帰係数 β の最良線形不偏推定量（Best Linear Unbiased
Estimator，BLUE）β̂ = (M  C −1 M )−1 M  C −1 Y であり，γ = m(s∗ ) − M  C −1 c∗ ，c0 = Var((s∗ ))とする（Cressie, 1993）．
BLUP は地球統計学の分野ではクリギング（Kriging）とも呼ばれている．クリギングは説明
変数を 3 種類に分け，p = 0, m(s) ≡ 0（期待値が既知の場合と同等）を単純クリギング（Simple，p = 1, m(x) = 定数 (未知) を通常クリギング（Ordinary Kriging），m(x) が一般の説明
Kriging）（Universal Kriging）と言う．
変数，β が未知の場合を一般クリギング
以下では単純クリギング Y (s) = (s) を考える．したがって（3.1）と（3.2）は各々Ŷ (s∗ ) = (c∗ ) C −1 Y

(3.3)
(3.4)

E[Y (s∗ ) − Ŷ (s∗ )]2 = c0 − (c∗ ) C −1 c∗になる．
次に C0 (h) を真の自己共分散関数，C1 (h) を我々が予測量を構成する際に想定する自己共分
散関数とする．また fi (λ)(i = 0, 1) を各々のスペクトル密度関数とする．真の BLUP と想定された自己共分散関数に基づく BLUP は（3.3）よりŶ (x∗ , Ci ) = (c∗i ) Ci−1 Y (i = 0, 1)となる．
サンプル数 n が増加するとき二つの予測量に関して
(3.5)

lim

EC0 [Y (s∗ ) − Ŷ (s∗ , C1 )]2=1∗∗
2
C0 [Y (s ) − Ŷ (s , C0 )]n→∞ E時空間大規模データに対する統計的解析法

61が成立する条件を考えたい． EC0 は真の自己共分散関数の下での期待値を表す．予測量の定義から任意の n に対して
（3.5）の比は 1 以上であるが，サンプル数が増加するに連れて想定された自己共分散関数に基づく BLUP の予測誤差が真の BLUP のそれと差がなくなることを等号は意味している．
以下では {Y (s)} は正規確率場に従うと仮定する．ただし予測誤差（3.4）は確率場の自己共分散関数のみで決まり，分布に関する他の条件は一切使用しないので，この仮定は一般性を失わない．もし Ci (h)(i = 0, 1) から導入される正規確率場の分布 Pi (i = 0, 1) が互いに「近ければ」
，関
係式
（3.5）
の成立が期待できる．確率分布の「近さ」を測るひとつの概念が同値性（equivalence）
である．
定義 1. Ω を確率空間，F を Ω 上で定義された σ-ﬁeld，Pi (= 0, 1) を確率分布とする．
P0 (A) = 0 ⇐⇒ P1 (A), A ∈ F .

が成立するとき，P0 と P1 は同値であると言う．
このとき以下の定理が成立する．
（3.5）が成立する．
定理 1.（Stein, 1990）確率分布 P0 と P1 が同値なとき，定理 1 は非定常確率場においても成立する．定常確率場に対しては確率分布の同値性が成立するための十分条件がスペクトル密度関数を用いて表現でき，以下の定理を得る．
定理 2.（Stein, 1993）以下の条件を仮定する．
(3.6)
(3.7)

0 < liminf||λ||→∞f0 (λ)/|φ(λ)|2 ≤ lim sup f0 (λ)/|φ(λ)|2 < ∞||λ||→∞

lim f1 (λ)/f0 (λ) = a(> 0)||λ||→∞

ここで φ(λ) は band-limited 関数である．band-limited 関数とはそのフーリエ変換のサポートがコンパクト集合となる関数を言う．このとき地点 s∗ が {si |i = 1, 2, . . . } の集積点であれば（3.5）
が成立する．
ここで定理 2 について補足しておく．
（3.6）
はあくまで技術的な仮定であり，これなしに定理2 が成立するか否かは未解決な問題である（Stein, 1999）
．また
（3.7）
は，値を予測したい地点 s∗が {si |i = 1, 2, . . . } の集積点のとき，予測量の精度は想定されたスペクトル密度関数の低周波の挙動には無関係で，高周波の挙動のみに依存することを意味している．
ここまでの議論においてはサンプル数 n に関わらず我々が想定する自己共分散関数 C1 (h) は固定されていた．実際のデータ解析においてはサンプル数に応じて C1 (h) を変化させていく方が自然である．この場合確率分布の近さを意味する概念として，近接性（contiguity）がある．
定義 2. {P0,n }，{P1,n } を確率空間の列 (Ωn , An ) において定義された確率分布の列とする．
任意の An ∈ An に対して，P0,n (An ) → 0 のとき P1,n (An ) → 0 が成立するならば，{P1,n } は{P0,n } に対して近接していると言う．また逆も同時に成立するとき， {P0,n } と {P1,n } は互いに近接していると言う．
このとき定理 2 は以下のように拡張できる．
定理 3.（Putter and Young, 2001） n → ∞ のとき EC0 [Y (s∗ ) − Ŷ (x∗ , C0 )]2 → 0 が成立し，また確率分布列 {P0 } と {P1,n } は互いに近接しているとする．

62

統計数理

第 60 巻

第 1 号 2012

このとき C1,n (h) を P1,n によって規定される自己共分散関数とすれば，P1 を P1,n に置き換えても定理 1 は成立する．
P0 および C0 (h) は各々真の確率分布，自己共分散関数であるから n には依存しないことを注意しておく．
定常確率場の場合には，近接性をスペクトル密度関数で表現すると以下のようになる．

定理 4.（Putter and Young, 2001） f1,n (λ) を自己共分散関数 C1,n (h) のスペクトル密度関数とする．
このとき {P1,n } が {P0 } に対して近接するための十分条件は
2

f1,n (λ) − f0 (λ)
lim sup n sup<∞
f0 (λ)n→∞λf1,n (λ)<∞
lim inf inf
n→∞ λ
f0 (λ)である．
逆に {P0 } が {P1,n } に対して近接するための十分条件は上式において f0 (λ) と f1,n (λ) を入れ替えた条件になる．
Covariance Tapering は真の自己共分散関数をサポートがコンパクト集合になる自己共分散関数によって代替し統計的解析を行う手法である．上述の予測問題に対する応用方法は以下のようになる．
いま Cθ (h) をサポートがコンパクト集合（Cθ(h) = 0, ||h|| > θ）
である自己相関関数 (Cθ (0) = 1)とする．真の自己共分散関数 C0 (h) に Cθ (h) をかけて
Ctap (h) = C0 (h)Cθ (h)を定義する．図 1 に例を示す．ここでは C0 (h)， Cθ (h) ともに等方型を仮定し，θ = 0.8 としている．Ctap (h) が自己共分散関数のときの共分散行列を Σtap = (σij,tap )(n × n) とすれば，
Ctap (si − sj ), ||si − sj || ≤ θ
σij,tap =0, ||si − sj || > θ

図 1.

Covariance Tapering.

時空間大規模データに対する統計的解析法

63となる．したがって Σtap の成分の多くは 0 となる “sparse” な行列であり，BLUP を構成する
際に必要となる逆行列にたいしては高速な計算手法が数値解析の分野において開発されている（Davis, 2006）．
いま {Y (s)} は C0 (h) が Matérn 族に属する定常確率場とし，我々が予測に用いる自己共分散関数は C1 (h) = Ctap (h) とする．
このとき以下の定理を得る．
定理 5.（Furrer et al., 2006） Cθ (h) を等方型の自己相関関数 Cθ (h) = Cθ (||h||) とし，fθ (λ)
をそのスペクトル密度関数とする．またある (≥ 0) と θ のみに依存する定数 M (θ) < ∞ が存在しM (θ)
0 < fθ (λ) ≤(1 + ||λ||2 )ν+d/2+
をみたすと仮定する．このとき（3.5）が成立する．
Matérn 族に属する自己共分散関数は自動的に定理 2 の仮定（3.6）をみたすことが知られている（Stein, 1999）．したがって（3.7）が成立することを示せばよい．
4.

最尤推定法に対する近似方法

大規模かつ不規則間隔で観測されるデータに対して最尤推定量を適用するのは，等間隔に観測されるデータにも増して共分散行列の逆行列や行列式の計算が複雑になり実用的とは言えな
い．より実用的で高速計算可能な近似方法が時空間領域および周波数領域において提案されている．本章ではこれらの方法のいくつかを紹介し，その優劣について議論する．
まず時空間領域において Vecchia（1988）
によって提案され，後に Stein et al.（2004）
によって一般化された方法を説明する．観測値 Y (si )(i = 1, . . . , n) が与えられたとき，Y n = (Y (s1 ), . . . , Y (sn ))
と置く．{Y (s)} が期待値 0 の正規定常確率場にしたがう場合，定数項を除去した θ の対数尤度関数は11
Ln (θ; Y n ) = − log det(Σn ) − Y n Σ−1
n Yn
22になる．ここで Σn は n × n 共分散行列で，その (i, j) 成分は σ(si , sj ; θ) = Cov(Y (si ), Y (sj )) とし，θ はパラメータ・ベクトルとする．
上述のように行列式 det(Σn ) や逆行列 Σ−1n の計算は複雑かつ time-consuming になる．そこでまず Y n を b 個の部分ベクトル Y jn (j = 1, . . . , b) に分割し，Y n = (Y 1n , . . . , Y bn ) と表す．また Y (jn) = (Y 1n , . . . , Y jn ) と置く．このとき Ln (θ; Y n ) は条件付き対数尤度関数の和Ln (θ; Y n ) = Ln (θ; Y 1n )bLn (Y jn |θ; Y (j−1,n) )j=2

によって表現できる．ここで Ln (Y jn |θ; Y (j−1,n) ) は Y (j−1,n) が与えられたときの Y jn の条件付き対数尤度関数である．
さらに S (jn) を Y (jn) の部分ベクトルとして，
Ln (θ; Y n ) ≈ Ln (θ; Y 1n )bLn (Y jn |θ; S (j−1,n) )j=2により対数尤度関数を近似し， Ln (θ; Y 1n ) bj=2 Ln (Y jn |θ; S (j−1,n) ) を最大にする θ を推定す
る．真の最尤推定量より計算は簡便になるがその近似精度は Y n の分割方法，S (jn) の選択などに依存する．

統計数理

64

第 60 巻

第 1 号 2012

周波数領域における推定方法として 2 つ紹介する．ひとつは Fuentes（2007）により提案された方法である．簡単のため d = 2 とする．まず t = (t1 , t2 ) を中心とする一辺の長さが ∆ の正方
形上で Y (s) を積分して，{X(t)} を
 
(4.1)
X(t) = ∆−2
Y (s)ds|ti −si |≤∆/2(i=1,2)

によって定義する．このとき {X(t)} も定常確率場になり，そのスペクトル密度関数はfX (λ) = ∆−4 |Γ(λ)|2 fY (λ)

である．ここで Γ(λ) = 2i=1 2 sin(∆λi /2)/λi とする．
次に {X(t)} を離散地点 ∆z(z ∈ Z 2 ) においてサンプリングした定常確率場 {X(∆z)} を考える．{X(∆z)} のスペクトル密度関数は
(4.2)
f∆X (λ) =fX (λ + 2πQ/∆), λ ∈ [−π/∆, π/∆]2 .
2Q∈Zとなる．{X(∆z)} に対して時系列解析における最尤推定法を近似する方法として用いられる
Whittle 推定法を適用する．いま ∆X m = (X(∆z 1 ), . . . , X(∆z m )) (z i ∈ Z 2 ) と置き，目的関数をIm (λ)L∗m (θ; ∆X m ) =dλlog f∆X (λ; θ) +f∆X (λ; θ)[−π/∆,π/∆]2i∆z j λ 2によって定義する． Im (λ) はピリオドグラム Im (λ) = | m| である．L∗m (θ; ∆j= X(∆z j )e
X m ) を最小にする θ をパラメータとして採用する．
ただし実際の計算では（4.2）の無限級数は適当な有限級数で近似する．また X(∆z) は観測できないので，1
Y (si )
(4.3)
X̃(∆z) =nz s ∈Di∆z

によって代替する．ここで D∆z は ∆z を中心とする一辺の長さ ∆ の正方形，nz はそれに含まれる si の個数とする．
もう一つは Matsuda and Yajima（2009）によって提案された方法である．ここでは以下のよ
うなサンプリング・スキームを考える．観測地点は
si = (A1 ui,1 , . . . , Ad ui,d ) ,i = 1, . . . , nとする．ここでui = (ui,1 , . . . , ui,d ) ,i = 1, . . . , n ,dは [0, 1] をサポートに持つ確率密度関数 g(x) から生成される i.i.d. サンプルとする．各辺をAj = Aj (k)(j = 1, . . . , )，サンプル数を n = nk とおき k → ∞ のとき発散すると仮定する．
このとき観測値のフーリエ変換およびピリオドグラムをd1

Jk (λ) = (2π)− 2 G−1/2 |Sk | 2 n−1knkYsp exp(−isp λ)

p=1

Ik (λ) = |Jk (λ)|2によって定義する．ここで Sk = [0, A1 ] × · · · × [0, Ad ]，|Sk | はその体積 |Sk | = A1 × · · · × Ad まd
た G = [0,1]d g(x)2 dx とする．Jk (λ) は (2π)− 2 |Sk |−1/2 S Ys exp(−is λ)ds を近似する統計量kである．ただし G は未知なのでカーネル関数などを用いてノンパラメトリックに推定する．

時空間大規模データに対する統計的解析法

65

目的関数は Whittle 推定量と同じアイディアに基づきIk (λ)
dλlog{f (λ : θ) + ck (θ)} +(4.4)
Lk (θ) =f(λ:θ) + ck (θ)
D
によって定義する．ここで D は λ ∈ D → −λ ∈ D をみたす対称なコンパクト集合とする．ck (θ)
はピリオドグラムのバイアスを考慮して必要となる θ の関数である．Lk (θ) を最小にする推定量は k → ∞ のとき一致性・漸近正規性をみたす．
ここで 3 つの方法を比較したシミュレーション結果を紹介する（Matsuda and Yajima, 2009）．
自己共分散関数は
 1/2 ν 1/2 
2ν r
2ν r
σ1
C(h1 , h2 ) = σ0 I(h1 = 0, h2 = 0) + ν−1
Kν
2
Γ(ν)
ρ
ρ
r 2 = [β{h1 cos(θ) − h2 sin(θ)}]2 + [β −1 {h1 sin(θ) + h2 cos(θ)}]2 .

によって定義する．右辺第 1 項の σ0 I(h1 = 0, h2 = 0) はナゲット効果（nugget eﬀect）と呼ばれ，
hi (i = 1.2) → 0 のとき，自己共分散関数がジャンプする場合を考慮した量である．第 2 項はパラメタライゼーションは異なるが，2 章で説明した Matérn 族に属する自己共分散関数である．
一方 r は単純な距離とはせず β, θ は非等方性を表現するパラメータである．まず座標 (h1 , h2 )を角度 θ だけ回転させl1
cos θ − sin θ
h1
=
sin θ
cos θ
l2
h2
を新たな座標とする．次に座標 (l1 , l2 ) において l1 は β だけ拡大し，l2 は逆に β だけ縮小する．
その結果距離の等高線は楕円
r 2 = β 2 l12 + β −2 l22 .
になる．β = 1, θ = 0 のとき γ(h1 , h2 ) は等方型モデルになる．
このときスペクトル密度関数はf (λ1 , λ2 ) = φ(α2 + ω 2 )−(ν+1)となる．ここでω 2 = [β −1 {λ1 cos(θ) − λ2 sin(θ)}]2 + [β{λ1 sin(θ) + λ2 cos(θ)}]2√
α = 2 ν/ρ, φ = να2ν σ1 /πとする．真のパラメータは
ν = 3.0, σ1 = 10, ρ = 3.0, σ0 = 3.0, β = 0.85, θ = 0.3と置く．サンプル地点は [0, 40] × [0, 40] 上の一様分布にしたがって抽出し，サンプル数は 2000，（4.4）の積分領域は D =繰り返し回数は 200 回とした．Matsuda and Yajima（2009）の目的関数
2[−3π, 3π] とした．
以下の表がその結果である．表 1 の上段が平均値，下段の括弧内の数値が標準偏差である．
最後の 2 行は 1 回当たりの計算に要する CPU タイムの平均値とその標準偏差である．表 2 は
等方性 (β = 1, θ = 0) を帰無仮説とする仮説検定に尤度比検定を適用した結果である． χ2 分布に基づく名目的な有意水準を 0.05 に取ったときの，実際の水準（size）と検出力（power）を示している．
表からは Stein et al.（2004）の方法が ν の推定を除いて最良であり，また尤度比検定 における検出力も一番高い．ただし Matsuda and Yajima（2009）の方法は ν の推定に関しては最良で

統計数理

66

第 60 巻

第 1 号 2012

表 1.

推定量の平均値と標準偏差．

表 2.

尤度比検定の水準と検出力．

ある．理由は自己共分散関数よりスペクトル密度関数の方が ν の値が異なるとその形状の相違がより顕著になることにあると思われる．Fuentes（2007）の方法は推定・検定とも Stein et（4.2）の無限級数を有限の級
al.（2004）あるいは Matsuda and Yajima（2009）ほど良くはない．
数により，さらに
（4.1）
の積分を
（4.3）
の和によって各々近似することに起因すると考えられる．
Matsuda and Yajima（2009）の優れた点は計算速度の速さである．ちなみに平均計算時間（8.7秒）は Stein et al.（2004）の方法（189.2 秒）の 5％しか要しない．
5.

潜在過程モデル

時空間データはしばしば観測誤差を伴い，直接には観測できない潜在的な確率場の推測が目的になる場合がある．本章では潜在的な確率場の推測を高速に行うために提案された 2 つの方法を紹介する．まず観測される確率場 Z(s) を(5.1)

Z(s) = Y (s) + U (s)によって定義する．ここでは潜在する確率場 {Y (s)} が 3 章で述べた回帰モデル(5.2)

Y (s) = m(s) β + (s)にしたがう．ただし (s) が非定常確率場になることも考慮してその自己共分散関数をCov((s), (s )) = C(s, s )と置く．{U (s)} は観測誤差などを意味する期待値 0 分散 τ 2 v(s) の互いに無相関な確率変数列とし，v(s) は既知の関数とする．

時空間大規模データに対する統計的解析法

67

観測値ベクトルを Z = (Z(s1 ), . . . , Z(sn )) と置く．C, V を各々 C = (C(si , sj ))，V = diag(v(s1 ), . . . , v(sn )) によって定義される n × n 行列とすれば，Z の共分散行列は
Σ = C + τ 2V

(5.3)

になる．
データ解析を高速化するために，多変量解析における因子分析と同様の発想に基づき共分散行
列 C を次元を縮約した行列により表現する方法がある．ひとつは Banerjee et al.（2008）
によって提案された正規予測過程（Gaussian Predictive Process）
であり，他方は Cressie and Johannesson（2008）によって提案された固定階数クリギング（Fixed Rank Kriging）である．
正規予測過程はデータ観測において主要と思われる地点を m(< n) 個まず選び，D∗ = (s∗1 , . . . ,∗(s)} は (s∗i )(i = 1, . . . , m) が与えられたときの BLUPsm ) とする．予測過程 {˜˜(s) = c(s) (C ∗ )−1 ∗によって定義する．ここで C ∗ は m × m 行列，c(s)，∗ は各々 m 次元ベクトルでC ∗ = (C(s∗i , s∗j ))c(s) = (c1 , . . . , cn ) , ci = Cov((s), (s∗i ))∗ = ((s∗1 ), . . . , (s∗m ))

とする．（5.2）において (s) を ˜(s) で置き換えてZ(s) = m(s) β + ˜(s) + U (s)(5.4)

を正規予測過程と呼ぶ．（5.1）と（5.4）は同じ Z(s) により表現しているが，異なった確率場である．
一方固定階数クリギングでは（5.2）において (s) = B  (s)η と置く．ここで η は m(< n) 次元
確率ベクトルでその共分散行列を G とする．B(s) = (B1 (s), . . . , Bm (s)) は m 個の基底関数の集合（必ずしも直交しない）である．このとき（5.3）はΣ = BGB  + τ 2 V(5.5)となる．ここで B は n × m 行列でその (i, j) 成分は B j (si ) とする．B およびGの選択により，多様な非定常確率場の共分散行列を表現できる．また（5.5）の Σ の逆行列は
Σ−1 = (τ 2 V )−1 − (τ 2 V )−1 B(G−1 + B  (τ 2 V )−1 B)−1 B  (τ 2 V )−1である．V は対角行列なので，実質的に Σ−1 の計算は m × m 行列 G−1 + B  (τ 2 V )−1 B の逆行列の計算に帰着する．
6.

今後の課題

最後に本章では今後解決すべき課題について列挙する．
（i）Covariance Tapering については，まず実際のデータ解析では真の自己共分散関数は未知であり，推定量を代入した場合の影響を明らかにしなければならない．
次に Taper のレンジ θ をサンプル数に依存させてどの様に取るべきかについて実用上の指針が必要である．次に Matérn 族に代表される等方型以外の非等方型あるいは分離型自己共分散関数に対してその理論的性質を導く必要がある．たとえば分離型モデルにおいては，自己共分散関数を 2 つに分割した場合，C(h) から生成される観測値間の共分散行列は Ci (h̃i )(i = 1, 2)から生成される共分散行列のクロネッカー積によって表現できる．したがってその逆行列も

統計数理

68

第 60 巻

第 1 号 2012

各々の逆行列のクロネッカー積になり，大規模データの処理に都合がよいモデルである．いまC1 (h̃1 ) を空間自己共分散関数， C2 (h̃2 ) を時間自己相関関数とすれば，空間自己共分散関数を
tapering すれば，さらに高速計算が可能になる．この場合に定理 4，5 が一般化できるかまだ未解決である．
最後に対数正規確率場など正規定常確率場の非線形変換により得られる定常確率場への拡張が考えられる．いま {Y (s)} は期待値 0，分散 1，自己共分散関数は Matérn 族に属する正規確率場とする．実数値関数 T (x) : R → R を用いて，新たな確率場 {Z(s)} をZ(s) = T (Y (s)) ,によって定義する．たとえば T (x) = exp(x) とおけば，{Z(s)} は対数正規確率場になる．φ(x)2を標準正規分布の密度関数 φ(x) = √12π e−x /2 としたとき，T (x) が ∞T (x)2 φ(x)dx < ∞ ,−∞をみたすならば，Z(s) はエルミート多項式の無限和
Z(s) =∞
aj Hj (Y (s)) ,j=1

によって表現できる．ここで Hj (x)(j = 1, 2, . . . ) は j 次のエルミート多項式であり，右辺は平均二乗収束による無限和である．Z(s) の自己共分散関数とスペクトル密度関数はCZ (||h||) =∞
a2j j!CY (||h||)j
j=1

∞

(j)
fz (||λ||) =a2j j!fY (||λ||) ,j=1となる．ここで fY(j) (||λ||) は fY (||λ||) の j 回結合関数である．
このような確率場についても定理 5 が成立すると予想される（Hirano, 2011）．
（ii）本稿で紹介した 3 つの尤度関数に対する近似法についても各々今後考察すべき課題がある．まず Stein et al.（2004）に関しては，前述のごとくシミュレーションにおいては良好なパ
フォーマンスを示すが，その理論的正当性は導出されていない．また観測値の順序づけ，グルー
プへの分割方法及びその個数について実際的なガイドラインが必要である．
Fuentes（2007）については正方形の一辺の長さ ∆ の決定方法および（4.2）の無限級数を有限級数で切断する方法についてより詳細に議論する必要がある．
最後に Matsuda and Yajima（2009）においては，推定量の精度は積分領域 D の選択に依存する．サンプルを無駄にしないためにはサンプル数が増大するとともに領域 D も拡張していく
ことが自然に思われるが，どのような速度で増大させれば一致性および漸近正規性が依然として成立するかまだ未解決である．
（iii）
正規予測過程については地点 D∗ = (s∗1 , . . . , s∗m ) およびその個数 m の選択方法の開発が必要となる，一方固定階数クリギングについては基底関数 B(s) = (B1 (s), . . . , Bm (s)) の選択方法とその個数 m の選択方法の開発が必要となる．同時に両方法の長所，短所，互いの優劣なども今後解明されなければならない．
謝辞

初稿を改訂するにあたり，査読者から大変有益なご助言をいただいた．ここに謝意を表したい．

時空間大規模データに対する統計的解析法

69

参 考 文 献Banerjee, S., Gelfand, A. E., Finley, A. O. and Sang, H.（2008）.

Gaussian predictive process models

for large spatial data sets, Journal of the Royal Statistical Society, B70, 825–848.
. Statistics for Spatial Data, rev. ed., Wiley, New York.
Cressie, N.（1993）
Cressie, N. and Johannesson, G.（2008）. Fixed rank kriging for very large data sets, Journal of theRoyal Statistical Society, B70, 209–226.
. Direct Methods for Linear Sparse Systems, SIAM, Philadelphia.
Davis, T. A.（2006）
Du, J., Zhang, H. and Mandrekar, V. S.（2009）. Fixed domain asymptotic properties of tapered
maximum likelihood estimators, The Annals of Statistics, 37, 3330–3361.
. Approximate likelihood for large irregularly spaced spatial data, Journal of the
Fuentes, M.（2007）
American Statistical Association, 102, 321–331.
. Covariance tapering for interpolation of large spatial
Furrer, R., Genton, M. and Nychka, D.（2006）
data sets, Journal of Computational and Graphical Statistics, 15, 502–523.
Guttorp, P. and Gneiting, T.（2006）. Studies in the history of probability and statistics XLIX, On
the Matérn correlation family, Biometrika, 93, 989–995.
. Random Fields on a Network, Modeling, Statistics and Applications, Springer,
Guyon, X.（1995）
New York.
. Covariance tapering for prediction of large spatial data sets in transformed random
Hirano, T.（2011）
ﬁelds（in preparation）.
. Covariance estimation for likelihoodKaufman, C. G., Schervish, M. J. and Nychka, D. W.（2008）
based estimation in large spatial data sets, Journal of the American Statistical Association,
103, 1545–1555.
Matérn, B.（1947）. Metoder art Uppskatta Noggranhetten vid Linje- och Provytetaxering, Stockholm:
Medd Staten Skogsforskningsinstitut, 36, no. 1，（In Swedish with substantial English summary）．
Matérn, B.（1960）.

Spatial Variation-stochastic Models and Their Application to Some Problems

in Forest Surveys and Other Sampling Investigations, Stockholm: Medd. Statens Skogsforskningsinstitut, 49, no. 5.
Matérn, B.（1986）. Spatial Variation, 2nd ed., Lecture Notes in Statistics, Springer, Berlin.
Matsuda, Y. and Yajima, Y.（2009）. Fourier analysis of irregularly spaced data on Rd , Journal of
the Royal Statistical Society, B71, 191–217.
. On the eﬀect of covariance function estimation on the accuracy
Putter, H. and Young, G. A.（2001）
of kriging predictors, Bernoulli, 7, 421–438.
. Uniform asymptotic optimality of linear prediction of a random ﬁeld using an
Stein, M. L.（1990）
incorrect second-order structure, The Annals of Statistics, 18, 850–872.
. A simple condition for asymptotic optimality of linear prediction of random
Stein, M. L.（1993）
ﬁelds, Statistics and Probability Letters, 17, 399–404.
Stein, M. L.（1999）. Interpolation of Spatial Data: Some Theory for Kriging, Springer, New York.
. Approximate likelihood for large spatial datasets,
Stein, M. L., Chi, Z. and Welty, L. J.（2004）
Journal of the Royal Statistical Society, B66, 275–296.
Su, Y., Li, B. and Genton, M.（2011）. Geostatistics for large data sets, Space-Time Processes andChallenges Related to Environmental Problems: Proceedings of the Spring School “Advances
and Challenges in Space-time Modelling of Natural Events”（eds. E. Porcu, J. M. Montero, Springer, Berlin（to appear）.
and M. Schlather）
Vecchia, A. V.（1988）. Estimation and model identiﬁcation for continuous spatial processes, Journalof the Royal Statistical Society, B50, 297–312.

70

統計数理

第 60 巻

第 1 号 2012

Wang, D. and Loh, W-L.（2011）. On ﬁxed-domain asymptotics and covariance tapering in Gaussian
random ﬁelds models, Electronic Journal of Statistics, 5, 238–269.
. Correlation Theory of Stationary and Related Random Functions, Vol. I, SprinYaglom, Y.（1987）
ger, New York.
矢島美寛 （2011）. 時系列解析から時空間統計解析への展望，日本統計学会誌，41, シリーズ J, 219–244.

Proceedings of the Institute of Statistical Mathematics Vol. 60, No. 1, 57–71 (2012)

71

Statistical Analysis of Large Spatio-temporal Data Sets
Yoshihiro Yajima and Toshihiro Hirano
Graduate School of Economics, The University of Tokyo

We review various time-saving spatio-temporal statistical methodologies and discussproblems to be solved in future. First we consider covariance tapering for the best linear
unbiased estimator (BLUP), which is called kriging in geostatistics. Secondly we consider
likelihood approximation in both spatio-temporal and frequency domains. Thirdly wedescribe latent process models which reduce the number of the parameters substantiallyso that it is able to analyze large spatio-temporal data sets within feasible computational
time. Finally we discuss open problems to be solved in future

Key words: Spatio-temporal statistical analysis, large spatio-temporal data sets, covariance tapering,latent processes.