特集「時空間統計解析：新たなる分野横断的展開」［原著論文］

統計数理（2012）
第 60 巻 第 1 号 131–148
c 2012 統計数理研究所
ブロック行列による正規確率場の逐次シミュレーション法について††

間瀬 茂 ・持田 信行 ・村山 靖洋†（受付 2011 年 6 月 29 日；改訂 2012 年 1 月 27 日；採択 1 月 27 日）

要旨

2 次定常正規確率場の理論的に正確な（少なくとも有界なレンジの共分散関数に対して）シ
ミュレーション法としては，共分散行列のコレスキ分解による方法が直接的であるが，共分散行列のサイズが大きくなると数値的困難が生じる．一方で，標準的シミュレーション方法である FFT（高速フーリエ変換）法を用いた方法は，規則的な格子点上のシミュレーションにしか使えない．この論文では，ブロック化された多変量正規変数の条件付き平均と条件付き分散の一般的な逐次公式を導き，それを用いて正規確率場をブロック毎に逐次的にシミュレーションする方法を提案する．この方法は少なくとも理論的には FFT 法よりも計算量が少ない．

キーワード： 正規確率場，逐次的条件付き平均・分散公式，ブロック行列，逐次シ
ミュレーション．
1.

序

地球統計学
（Chilès and Delﬁner, 1999; 間瀬, 2010; Wackernagel, 2003）を典型として，確率場は連続的空間で観測されるデータのモデルとして重要である．特に条件付きシミュレーション
法はデータ観測地点でデータ値と一致する確率場を多数シミュレーションすることにより空間予測を行う．Matheron による巧妙な方法である確率場の条件付きシミュレーション法（Chilès
et al., 1999, pp. 465–468）を用いれば，クリギング予測と確率場の無条件シミュレーションから条件付きシミュレーションを構成できる．幾つかのアドホックな手法を除けば，確率場のシミュレーションは与えられた共分散関数を持つ 2 次定常（結果として強定常）
正規確率場を指定された位置の集合上で生成する．
2 次定常確率場のシミュレーション法は，与えられた有限個の位置での確率場の値を同時生成する．そうした方法としては，（1）Matheron による古典的な turning band 法（Matheron,（2）共分散行列のコレスキ分解を用いた直接法
1973; Chilès and Delﬁner, 1999, pp. 472–477），（Davis, 1987; Chilès and Delﬁner, 1999, pp. 464–465），（3）そして離散フーリエ変換と FFT
アルゴリズムを用いた方法（Pardo-Iguzquiza and Chica-Olmo, 1993; Chilès and Delﬁner, 1999,，が有名である．
（1）は任意位置集合での確率場を生成できるが，本質的に近似共
pp. 496–504）
分散構造を用い，相対的に計算速度は遅い．また低次元空間に於ける正規確率場をシミュレー
ションできることが前提となる．
（2）も任意位置集合での確率場を生成できるが，共分散行列全体のコレスキ分解を用いるため，そのサイズと分解速度が制約になる．
（3）は相対的に高速（位†

東京工業大学大学院 情報理工学研究科：〒152–8550 東京都目黒区大岡山 2–12–1 W8–28

132

統計数理

第 60 巻

第 1 号 2012

置の数を n とすると計算量は O(n log n)）であるが，位置集合は規則的である必要があり，また範囲が無限の共分散関数に対してはそれを有限範囲な関数で近似する必要がある．
この論文では（2）の手法を拡張し，正規確率ベクトルの条件付き平均と分散の逐次的公式を用い，領域を分割し確率場を部分領域毎に逐次的にシミュレーションする方法を提案する．正規ランダムベクトルを条件付き平均と分散公式を用いて逐次的に生成するというアイデア自体は古典的であるが，この論文では部分領域に対応してブロック化された共分散行列を用い，各部分領域上の正規確率場の条件付き平均と分散を逐次的に更新する公式を与える．この手法では，原理的には位置集合は任意で良く，計算量は O(n) になる．一方で条件付き平均と分散を逐次的に計算する分，実際の計算速度はやはり部分領域とその中の位置集合のサイズに依存す
る．更に部分領域とその中の位置集合が合同な場合，条件付き共分散行列が極限を持つことを示し，極限条件付き共分散行列を予め計算し，それを用いてシミュレーション速度を早める可能性と，その理論的正当性についても検討する．最後に提案手法を用いたシミュレーションの
数値実験結果を紹介する．
1.1 2 次定常確率場
d 次元空間を径数集合に持つ確率変数の集合 Z = {Z(x); x ∈ Rd } を確率場と呼ぶ．統計モデ
ルとしての確率場には以下の 2 次定常性の仮定が標準である：（1）平均値は一定値 E{Z(x)} = µ，（2）共分散は位置の差分 h のみに依存 C(h) = C{Z(x + h), Z(x)} する．
位置の差分 h のみに依存する関数 C(h) を確率場の共分散関数と呼び，次のような基本的性質を持つ：（1）偶関数： C(−h) = C(h)，（2）有界性：|C(h)| ≤ C(0) = V {Z(x)}，（3）非負値定値性： すべての c1 , c2 , . . . , cn で ni,j=1 ci cj C(xi − xj ) ≥ 0．
以下で用いる共分散関数は狭義正定値であるもののみを用いる．共分散関数 C(h) がノルム||h|| のみに依存 C(h) = C(||h||) するとき確率場は等方的であるといわれる．ある r0 > 0 があり，r ≥ r0 ならば C(r) = 0 となる時，共分散関数は有限なレンジを持つといい，そうした r0の最小値を C(·) のレンジと呼ぶ．地球統計学で用いられる代表的な等方的共分散関数モデルには以下のようなものがある：（1）指数型モデル：Cexp (r) = b exp(−r/a)（a, b > 0）．
（2）球型モデル：有限なレンジを持つ．
b(1 − 3r/(2a) + r 3 /(2a3 )) (r ≤ a),
Csph (r) =（その他）．
0（3） Gaussian モデル：Cgau (r) = b exp(−r 2 /a2 )，（a, b > 0）．
（4） Gneiting モデル：有限なレンジを持つ．
b(1 + 8sr + 25(sr)2 + 32(sr)3 ) (0 ≤ r ≤ a/s),
Cgntg (r) =（その他）．
0
ここで s = 0.301187465825 . . . ．このモデルは d ≤ 3 でのみ使え，6 階微分可能である．
1.2 ブロック行列に関する結果
この節では，参考のために，以下の正規確率場の逐次シミュレーション法の議論で必要にな

ブロック行列による正規確率場の逐次シミュレーション法について

133

る，主としてブロック行列に関する結果を，必ずしも直接引用しないものを含めまとめておく．
T ，U ，V ，W をそれぞれサイズ m × m，m × n，n × m，n × n の行列とし，ブロック行列
 U
A = VT W
を考える．
補題 1.（ハーヴィル, 2007, 定理 13.3.8）もし T が非特異であれば
  W V −1 
 = |T ||W − V T −1 U | .
|A| = 
U T
補題 2.（ハーヴィル, 2007, 定理 8.5.5）ブロック下三角行列 B =非特異行列であるときに限って非特異であり B −1 =T O T −1
O−W −1 V T −1 W −1

V W

は T と W が共に．

補題 3.（ハーヴィル, 2007, 定理 8.5.10，定理 8.5.11） T が非特異であれば A は Q = W −V T −1 U（T のシューアの補元）が非特異なときに限って非特異であり，その逆行列は A−1 = T −1 +T −1 U Q−1 V T −1 −T −1 U Q ．
Q−1−Q−1 V T −1
補題 4. もし A が対称ならば T のシューアの補元 Q = W − V T −1 U も対称である．
証明．A が対称ならば T t = T ，U t = V ，V t = U ，W t = W である．従って
Qt = W t − (V T −1 U )t = W − U t (T t )−1 V t = W − V T −1 U = Q .

補題 5.（ハーヴィル, 2007, 定理 14.8.4） T は正定値で非特異行列と仮定する．この時 A が正定値ならば T のシューアの補元も正定値である．
補題 6.（ハーヴィル, 2007, 系 14.2.12）正定値行列の任意の主部分行列は正定値である．
補題 7.（ハーヴィル, 2007, 定理 14.9.5） A = {aij } を n × n 対称行列とし，k = 1, 2, . . . , n に対して⎤⎡
a11 a12 · · · a1k⎥⎢a⎢ 21 a22 · · · a2k ⎥Ak = ⎢
..
.. ⎥
..
⎥
⎢ ..
.
⎣ .
.
. ⎦
ak1 ak2 · · · akk
を次数 k の A の首座部分行列とする． A のすべての首座部分行列 A1 , A2 , . . . , An の行列式が全て正のとき，かつそのときに限って A は正定値である．
補題 8.（ハーヴィル, 2007, 補助定理 14.9.6） An−1 を n × n 対称行列 A = {aij } の n − 1 次
首座部分行列とする．A が正定値である必要十分条件は An−1 が正定値かつ |A| > 0 であることである．
補題 9.（コレスキ分解，ハーヴィル, 2007, 定理 14.5.11）任意の n × n 対称正定値行列 A に対して，A = LLt を満たす正の対角要素を持つ可逆な下三角行列 L が一意に存在する．
以下では対称正定値行列 A のコレスキ分解の下三角行列を Ch(A) と表すことにする．多変量正規分布の条件付き平均と分散については，次の良く知られた結果が基本となる．
補題 10. n + m 次元変量正規確率ベクトル X = (X t1 , X t2 )t の平均ベクトルと共分散行列を

統計数理

134

第 60 巻

第 1 号 2012

それぞれ
µ = (E{X 1 }t , E{X 2 }t )t = (µt1 , µt2 )t ,
 

Σ11
C{X 1 , X 1 } C{X 1 , X 2 }=Σ=
C{X 2 , X 1 } C{X 2 , X 2 }
Σ21

Σ12
Σ22



と置く．条件 X 1 = x1 の下での X 2 の条件付き平均ベクトルと共分散行列は，それぞれ次のようになる：
E{X 2 | x1 } = µ2 + Σ21 Σ−1
11 (x1 − µ1 ) ,C{X 2 | x1 } = Σ22 − Σ21 Σ−1
11 Σ12 .

特に条件付き分散は x1 に依存せず非ランダムな量になる．
2.

ブロック化された条件付き平均ベクトルと共分散行列

平均ゼロで共分散行列 Σ を持つ正規確率ベクトル X を生成する最も単純で直接的な方法
は，共分散行列のコレスキ分解 Σ ＝ LLt と標準正規分布に従う独立な確率ベクトル Y を用いて X = LY とすることである．実際 E{XX t } = LE{Y Y t }L = LLt = Σ となる．この節では共分散行列をブロックに分割し，逐次的に正規確率ベクトルをシミュレーションする方法を考える．
以下では正規確率ベクトルの（正定値）共分散行列を k × k ブロックに分割したものを
⎡

Σ11
⎢Σ
⎢ 21
Σk = ⎢
⎢ ..
⎣ .
Σk1

Σ12
Σ22
..
.
Σk2

···
···
..
.

⎤
Σ1k
Σ2k ⎥
⎥
.. ⎥
⎥
. ⎦

···

Σkk

とし

Σk+1 =

Σk
Bk+1

t
Bk+1

Σk+1,k+1



ここで Bk+1 = [ Σk+1,1···

Σk+1,2

Σk+1,k ]と置く．また Σk のコレスキ分解を Σk = Lk Ltk とし，下三角行列 Lk を Σk に準じてブロック化した各ブロックを Lij と置く．当然 Lk はブロック下三角行列（i < j ならば Lij = O）
である．
Lkk の逆行列を Kkk とする．また E と O はそれぞれ適当なサイズの単位行列と零行列とする（同じ式の中でもサイズは異なる可能性がある）．さらにいくつかの証明中では，その具体的な形が不要なブロック成分をシンボル ∗ で表現する．
命題 1. Lk は一貫性を持つ，つまり
LkO
Lk+1 =, ここで Yk+1 = [ Lk+1,1
Yk+1 Lk+1,k+1

Lk+1,2···

証明．Lk+1 の左上の n × n ブロック部分行列を L∗k とすると

 ∗

 ∗
t
Σk
Lk
Bk+1
O
Lk
Σk+1 == Lk+1 Ltk+1 =Bk+1 Σk+1,k+1
O
Yk+1 Lk+1,k+1

Lk+1,k ] .

t
Yk+1

Ltk+1,k+1

.

ブロック行列による正規確率場の逐次シミュレーション法について

R=

E O
O O

135

を用いると




O
E
= RΣk+1 Rt =
O
O

Σk
O

O
O



Σk
Bk+1



t
Bk+1

Σk+1,k+1

E
O

O
O



となる．従って RΣk+1 Rt = RLk+1 Ltk+1 Rt = (RLk+1 )(RLk+1 )t であり

 ∗
  ∗

E O
Lk
Lk O
O
RLk+1 ==O O Yk+1 Lk+1,k+1
O O
となるから



Σk
O


 ∗
O
Lk
t
t
= RΣk+1 R = (RLk+1 )(RLk+1 ) =
O
O

O
O



L∗k
O

O
O

t
.

従って Σk = L∗k (L∗k )t である．コレスキ分解の一意性（補題 9）から L∗k = Lk となるから，結局
 L
Lk+1 = Y k L Oとなり，一貫性が保たれる．
k+1

k+1,k+1

命題 2. Lk の i 行目のブロック Lki は以下の性質を持つ：


i−1

Lk1 = Σk1 (Lt11 )−1 , Lki = Σki −Lka Ltia (Ltii )−1(2 ≤ i ≤ k − 1) ,a=1

Lkk Ltkk = Σkk −k−1


Lka Ltka .

a=1

証明．簡単のために Bk = [ Σk1

Σk =

Σk−1
Bk

Σk2

···

Σk,k−1 ]，Yk = [ Lk1


t 

Lk−1
Bk
= Lk Ltk =
Σkk
Yk

O
Lkk



···

Lk2

Ltk−1
O

Lk,k−1 ] と置くと

t 

Yk
Ltkk

となり，i = 1 の時 Σk1 = Lk1 Lt11 ，Lk1 = Σk1 (Lt11 )−1 ，2 ≤ i ≤ k − 1 の時は i < a で Lia = O で
あるから
Σki =

k−1


Lka Ltia =

a=1

i−1


k−1


Lka Ltia + Lki Ltii +

a=1



i−1

Lka Ltia (Ltii )−1 .
Lki = Σki −a=i+1

Lka Ltia =i−1
Lka Ltia + Lki Ltii ,

a=1

a=1

そして i = k の時
Σkk = Yk Ykt + Lkk Ltkk =k−1
Lka Ltka + Lkk Ltkk ,a=1

Lkk Ltkk = Σkk −k−1


Lka Ltka .

a=1

命題 3. 任意の正整数 s を固定する．|k − j| > s である全ての k, j に対して Σkj = O であ
れば Lkj = O となる．
証明．k を一つ固定する．j = 1 の時 Σk1 = O であれば Lk1 = Σk1 (Lt11 )−1 = O となる．j = l の時 Σk1 = Σk2 = · · · = Σkl = O ならば Lk1 = Lk2 = · · · = Lkl = O が成り立つと仮定する．j = l + 1

統計数理

136

第 60 巻

第 1 号 2012の時 Σkj = O であれば
Lkj = (Σkj −l


Lka Ltja )(Ltjj )−1 = O(Ltjj )−1 = O

a=1となり，j = l + 1 のときも Lkj = O となる．従って数学的帰納法により証明が得られる．
有界範囲の共分散関数に対しては，命題 3 は次のように拡張できる．
命題 4. 任意の正整数 s を固定する．|k − j| > s である全ての k, j に対して Σkj = O であれば，Lkj は以下の関係を持つ： k < s + 2 の時


j−1

Lk1 = Σk1 (Lt11 )−1 , Lkj = Σkj −Lka Ltja (Ltjj )−1 (2 ≤ j ≤ k − 1) ,a=1

Lkk Ltkk = Σkk −k−1


Lka Ltka .

a=1

また s + 2 ≤ k の時は
Lkj = O

(1 ≤ j ≤ k − s − 1) ,Lk,k−s = Σk,k−s (Ltk−s,k−s )−1 ,


j−1

Lka Ltja (Ltjj )−1
Lkj = Σkj −(k − (s − 1) ≤ j ≤ k − 1) ,

a=k−s
k−1


Lkk Ltkk = Σkk −Lka Ltka .

a=k−s

証明．k < s + 2 の時 Σk に Σij = O となるブロックはないので，命題 2 から直接結論がでる．s + 2 ≤ k の時は以下の 3 つの場合に分けて証明される： 1 ≤ j ≤ k − s − 1 の時は命題 3 より明らか．j = k − s の時は 1 ≤ j ≤ k − s − 1 で Lkj = O に注意すると


k−s−1

Lk,k−s = Σk,k−s −Lka Ltk−s,a (Ltk−s,k−s )−1 = Σk,k−s (Ltk−s,k−s )−1 .
a=1

k − (s − 1) ≤ j ≤ k − 1 の時は 1 ≤ j ≤ k − s − 1 で Lkj = O に注意すると




j−1
j−1


Lkj = Σkj −Lka Ltja (Ltjj )−1 = Σkj −Lka Ltja (Ltjj )−1 .
a=1

a=k−s

そして j = k の時は 1 ≤ j ≤ k − s − 1 で Lkj = O より
Lkk Ltkk = Σkk −

k−1

a=1

Lka Ltka = Σkk −k−1


Lka Ltka .

a=k−s

命題 5. Lk は非特異である．
証明．Σk は非特異だから |Σk | = |Lk Ltk | = |Lk ||Ltk | = |Lk |2 > 0．
以下では Vk は Lk の逆行列とし， {Vij } をそのブロックとする． Vk にも Lk と同様の一貫性がある．

ブロック行列による正規確率場の逐次シミュレーション法について

命題 6.



Vk+1 =

Vk
Rk+1

O
Vk+1,k+1

,

ここで Rk+1 = [ Vk+1,1

Vk+1,2···137

Vk+1,k ] .

また
Vk+1,k+1 = L−1
k+1,k+1 ,Vk+1,j = −Vk+1,k+1kLk+1,a Vaj(1 ≤ j ≤ k) .

a=j

証明．Lk+1 =
L−1
k+1 =



 L

O
k
Yk+1 Lk+1,k+1

の逆行列は，補題 2 より

L−1
k
−i
−Lk+1,k+1 Yk+1 L−1
k



O
L−1
k+1,k+1となる．また 1 ≤ j ≤ k で Vk+1,j = −Vk+1,k+1


=

k

Vk

O

−L−1
k+1,k+1 Yk+1 Vk

L−1
k+1,k+1

= Vk+1

a=j Lk+1,a Vaj となることも同時にわかる．

有界範囲の共分散行列の場合は以下の命題が成り立つ．
命題 7. 任意の正整数 s を固定する．|k − j| > s である全ての k, j に対して Σkj = O であれば，{Vkj } は以下の関係を満たす：
Vkk = L−1
kk ,

k−1


Vkj = −Vkk

Lka Vaj(1 ≤ j ≤ k − 1) .

a=max{j,k−s}証明．k < s + 2 の時 Σk に Σij = O となるブロックは無く，任意の j で Lkj = O となるので Vkj は命題 6 から
Vkk = L−1
kk ,Vkj = −Vkk

k−1Lka Vaj .

a=js + 2 ≤ k の時 Σkj = O（1 ≤ j ≤ k − s − 1）であるから，命題 3 より Lkj = O（1 ≤ j ≤ k − s − 1）
となる．よって j = k の時 Vkk = L−1．
k−s ≤ j ≤ k − 1 の時 Vkj = −Vkk k−1
kka=j Lka Vaj ．そして1 ≤ j ≤ k − s − 1 の時
Vkj = −Vkk

k−1Lka Vaj = −Vkk

a=j

k−1


Lka Vaj .

a=k−s

以上の結果より証明が完了する．
3.

逐次シミュレーション法

以上の結果を用いて，ブロック行列を用いた正規確率場の逐次シミュレーション法を考える．
一般性を失うこと無く，考える確率場の平均は総てゼロと仮定してよい．考慮する領域 A が n個の互いに排反な副領域群 A1 , A2 , . . . , An に分割されているとする．各領域 Ai にはそれぞれpi 個の点群を考え，確率場のその各位置での値を与える確率ベクトル Xi を生成する逐次アルゴリズムを考える．次の二つの命題は，条件付き平均と分散を逐次的に計算する方法を与える．
命題 8. 条件付き平均 Mk は次の漸化式を持つ：
M1 = ゼロベクトル,

Mk =

k−1

a=1

Lka Kaa (Xa − Ma )(k ≥ 2) .

統計数理

138

第 60 巻

第 1 号 2012tttt証明．条件付き平均は Mk = Bk Σ−1（補題 10）．Σk Σ−1
k−1 [X1 X2 · · · Xk−1 ] と表されるk =t
Σk (Kk Kk ) = E であるから
 t



O
Kk−1 Rkt
Kk−1
Σk−1 Bkt
E=
t
Bk
Σkk
O
Kkk
Kkk
Rk


∗
∗
=
t
t
t
Kk−1 + (Bk Rkt + Σkk Kkk)Rk (Bk Rkt + Σkk Kkk)Kkk
Bk Kk−1tとなる．ここで Rk = [Kk1 Kk2 · · · Kk,k−1 ] と置いた．従って (Bk Rkt + Σkk Kkk)Kkk = E であり
−1
t
t
t
t
Kkk = Bk Rk + Σkk Kkk である．この結果を代入すると Bk Kk−1 Kk−1 + (Bk Rkt + Σkk Kkk
)Rk =
−1
−1
−1
−1
Bk Σk−1 + Kkk Rk = O となることから Bk Σk−1 = −Kkk Rk = −Lkk Rk となる．この結果を更ttttに式 Mk = Bk Σ−1
k−1 [X1 X2 · · · Xk−1 ] に代入し，命題 2 の二番目の関係式を用いれば，次の逐次公式を得ることが出来る：

k−1
 
k−1
k−1



−1
−1
Mk = −Kkk
Kkj Xj = −Kkk
Lka Kaj Xj
−Kkk
j=1

=

k−1
k−1

j=1

=

k−1


j=1

a=j


a−1

k−1


Lka Kaj Xj =Lka
Kaj Xj + Kaa Xa

a=j

a=1

j=1

Lka (−Kaa Ma + Kaa Xa ) = −

a=1

k−1


Lka Kaa (Ma − Xa ) .

a=1

命題 9. 条件付き分散 Vk は関係 Vk = Lkk Ltkk = Σkk −Ch(Vk ) である．
証明．Vk は Σk のブロック化 Σk =

Σ

k−1

Bk

t
Bk
Σkk



k−1

t
a=1 Lka Lka を満たす．特に Lkk =における Σk−1 のシューア補元 Σkk −t
Bk Σ−1（補題 4，5）．従って Σ−1=
k−1 Bk であり，非特異で対称正定値である
k
t
で Σ−1
k = Kk Kk を用いると


 t
∗
∗
Kk−1
t
=
K
Σ−1
=
K
=
k k
k
O
∗ Vk−1∗∗ Vk−1



．ここ


∗
.
t
Kkk
Kkk

t
t
よって Vk = (Vk−1 )−1 = (Kkk
Kkk )−1 = Lkk Ltkk ．更に命題 2 より Lkk Ltkk = Σkk − k−1
a=1 Lka Lkak−1tであるから Vk = Σkk − a=1 Lka Lka ．
Qtk
t
Kkk



∗

Kk−1
Qk

 
O
∗
=
Kkk
∗

以上をまとめると，確率ベクトル X1 , X2 , . . . , Xn , . . . の逐次的生成法を与える以下の定理を得る．
定理 1. W1 , W2 , . . . , Wn , . . . を（適当なサイズの）標準正規分布に従う，互いに独立な成分を持つ，相互にも独立な確率ベクトルとすると
X1 = L11 W1 ,Xk = Lkk Wk + Mk

(k ≥ 2) .

ここでM1 = ゼロベクトル,

Mk =

k−1

a=1

行列 Lki は次の関係から逐次的に定まる：
L11 = Ch(Σ11 ) ,Lka L−1
aa (Xa − Ma )(k ≥ 2) .

ブロック行列による正規確率場の逐次シミュレーション法について

Lk1 = Σk1 (Lt11 )−1 (2 ≤ k) ,


i−1

Lka Ltia (Ltii )−1
Lki = Σki −139(2 ≤ i ≤ k − 1) ,

a=1



k−1

Lka Ltka .
Lkk = Ch Σkk −a=1

4.

有限なレンジの場合

前節の結果は，共分散関数が無限のレンジを持つ場合も有効であるが，当然ながらそうした
場合，必要な計算量は繰り返し毎に増大する．この節では有限なレンジを持つ共分散関数の場
合に，前節の結果が簡略化されることを示す．
定義 1. |i − j| > s を満たす全ての i, j に対して相互共分散行列 Σij = O となるとき，正整
数 s をレンジと呼ぶ．
共分散関数が有限なレンジを持てばそうした s が存在するが，その具体的な値は領域とその中の点列の配置に依存する．例えば s = 1 であれば⎡⎤··· ······
O
Σ11 Σ12 O
⎢
⎥
..
..
⎢
⎥
.
⎢ Σ21 Σ22 Σ23
⎥
.
⎢
⎥
⎢
⎥
.
.
.
..
..
..
⎢O
⎥
Σ
Σ
32
33
⎢
⎥
⎢.
⎥
..
..
..
..
..
..
⎥
.
Σk = ⎢
.
.
.
.
.
⎢.
⎥
.
⎢
⎥
⎢ ..
⎥
.
.
.
.
⎢.
⎥
.
.
Σ
O
Σ
k−2,k−2
k−2,k−1
⎢
⎥
⎢.
⎥
.
⎢.
⎥
.
. Σk−1,k−2 Σk−1,k−1 Σk−1,k ⎦
⎣.
O
···
···
··· O
Σk,k−1Σkkのようにブロック化される．従ってこの場合次の結果を得る．
定理 2. 有限なレンジ s を持つとき，条件付き平均と分散は以下の漸化式を持つ：k ≥ 2 で
Mk =

k−1


Lka L−1
aa (Xa − Ma ) ,a=max{1,k−s}

Vk = Σkk −k−1


Lka Ltka .

a=max{1,k−s}

証明．Σkj = Lkj = O（1 ≤ j ≤ k − (s + 1)）であることから，命題 8，9 を用いて結論が得られる．
5.

極限条件付き分散を用いた正規確率場の逐次シミュレーション法

これまでに述べた正規確率場の逐次シミュレーション法は，一度に扱う行列のサイズを減少させているが，それでも特に条件付き分散を計算する手間が毎回必要になる．一方で領域とその中の点群が合同である場合，数値実験によればこの非ランダムな条件付き分散は，比較的速やかに収束することが観察される．一方，条件付き平均は本来ランダムな量であり収束性は見

統計数理

140

第 60 巻

第 1 号 2012られない．従って，まずこの極限条件付き分散を求め，それを最初から使うことが考えられる．
この節ではこうした極限条件付き分散を用いた正規確率場の逐次シミュレーション法とその理論的妥当性を考察する．この節では図 1 のように合同な領域列 A1 , A2 , . . . が並んでおり，各 Ai
中の点列も格子状に並んでいる必要はないが互いに合同であるとする．以下の命題 10 の証明の必要上，領域列と確率ベクトル列はそれぞれ . . . , A−2 , A−1 , A0 及び . . . , X−2 , X−1 , X0 と延長されているとする．
命題 10. 条件付き共分散行列の列 Vk は，確率 1 で（非ランダムな）極限 V∞ を持つ．
証明．まず，必ずしも正規型でない強定常 2 次確率場を考える． Xn の i 番目の成分を
Xn:i と置くと，条件付き共分散行列 C{Xn | X1 , X2 , . . . , Xn−1 } の (i, j) 成分は E{Xn:i Xn:j |X1 , X2 , . . . , Xn−1 } と
E{Xn:i | X1 , X2 , . . . , Xn−1 }E{Xn:j | X1 , X2 , . . . , Xn−1 }の差に分解できるので，全ての i, j で
E{Xn:i Xn:j | X1 , X2 , . . . , Xn−1 },E{Xn:i | X1 , X2 , . . . , Xn−1 }のそれぞれが確率 1 で収束することを示せばよい． E{Xn:i | X1 , X2 , . . . , Xn−1 } の場合の証明を与える．
E{X0:i | X−1 , X−2 , . . . , X−(n−1) }（n = 1, 2, . . . ）はマルチンゲールでありその L1 ノルムは有界
E|E{X0:i | X−1 , X−2 , . . . , X−(n−1) }| ≤ E{|X0:i |}である．従って Doob のマルチンゲール前向き収束命題（Williams, 1991, p. 112）より E{X0:i |X−1 , X−2 , . . . , X−(n−1) } は確率 1 で有限な極限を持つことがわかる．同様な議論から E{X0:i X0:j |X−1 , X−2 , . . . , X−(n−1) } も確率 1 で有限な極限を持つ．結局，条件付き共分散行列 C{X0 |X−1 , X−2 , . . . , X−(n−1) } は確率 1 で有限な極限行列を持つ．
他方で，2 次定常な正規確率場であれば同時に強定常性を持つから， C{X0 | X−1 , X−2 , . . . ,X−(n−1) } と C{Xn | X1 , X2 , . . . , Xn−1 } は同じ確率分布に従う．さらに条件付き分散は非ランダ
ム（補題 10）でもあるから，上で議論した確率 1 の収束は非ランダムな量に対する収束を意味す
る．以上から結論が得られる．
次に s = 1 の場合に，極限行列 V∞ を数値的に求める方法を考える．この時 Lk,k−2 = 0 であ−1
ることに注意すると Vk = Σkk − Lk,k−1 Ltk,k−1 = Σkk − Σk,k−1 Vk−1
Σtk,k−1 となり，V∞ が存在す

図 1.

合同な領域列と対応する確率ベクトル列（上）と合同な点配置（下）．

ブロック行列による正規確率場の逐次シミュレーション法について141ることを考慮すると−1 t
V∞ = Σ11 − Σ21 V∞
Σ21

(5.1)

となり，極限条件付き分散に関する方程式が求まる．式 (5.1) の形の行列方程式の解法としては以下のものが考えられる（Fital-Akelbek, 2009; Benner and Faßbender, 2007; Varga, 2009）：（1）定義通りに Vk の数値的極限を求める，（2） ﬁxed-point iteration，（3） Meini’s algorithm，（4） doubling algorithm，（5） butterﬂy sz algorithm．
以下では L∗2 = Σ31 (Lt∞∞ )−1 と置き，L∗1 を行列方程式 L∗1 = (Σ21 − L∗2 (L∗1 )t )(Lt∞∞ )−1 の解とする．ここで L∞∞ = Ch(V∞ ) である．
ﬁxed-point iteration による式 (5.1) の解法は
V0∗ = Σ11 ,∗Vk∗ = Σ11 + Σ21 (Vk−1)−1 Σt21(k ≥ 2)となる．しかしこの方法よりは
V0∗ = Σ11 ,∗Ak−1 = Ch(Vk−1),

Bk = Σ21 A−1
k−1 ,Vk∗ = Σ11 + Bk Bkt(k = 1, 2, . . . )とした方が良い場合が多いようである．
一方，Meini’s algorithm（Fital-Akelbek, 2009, pp. 64–83）
による式 (5.1) の解法はつぎのようになる：
V0∗ = Q0 = Σ11 , A0 = Σ21 ,∗t= Xn − An Q−1
Vn+1
n An ,tt−1
Qn+1 = Qn − An Q−1
n An − An Qn An ,An+1 = An Q−1
n An(n = 0, 1, . . . ) .

次に，Vk の代わりに最初から極限条件付き分散 V∞ を用いて確率場を生成した場合，共分散構造が保たれているかどうかを理論的に検討する．
s = 1 の場合 V1 = Σ11 , L11 = Ch(V1 ) であり，k ≥ 2 では
∗
Vki
=O

(1 ≤ i ≤ k − 2) ,L∗k,k−1 = Σ21 (Lt∞∞ )−1 ,となり，従って⎡L11⎢⎢ Σ21 (Lt∞∞ )−1⎢⎢OL∗k = ⎢
⎢
..
⎢
.
⎣
O

Vk∗ = V∞ ,OL∞∞
Σ21 (Lt∞∞ )−1..
.
OL∗kk = L∞∞

···
···
···
..
.

O
O
O
..
.

···

L∞∞

⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦

統計数理

142

第 60 巻

第 1 号 2012となる．これより，極限条件付き分散を用いた場合の領域全体の共分散 Σ∗k = L∗k (L∗k )t は
⎡
⎤
t
Σ11
L11 L−1···
OO∞∞ Σ21⎢⎥⎢ Σ21 (Lt∞∞ )−1 Lt11
⎥
Σ11
···
O
O
⎢
⎥
⎢
⎥
..
..
.
.
.
.
.
.
⎢
⎥
.
.
.
.
.
⎢
⎥
⎢
⎥
−1 t
t
Σ21
O
O
· · · Σ21 V∞ Σ21 + V∞
⎣
⎦
−1 t
O
O
···
Σ21
Σ21 V∞ Σ21 + V∞となる．この式に (5.1) を代入して⎡
Σ11⎢ Σ (Lt )−1 Lt⎢ 21 ∞∞
11⎢⎢O⎢Σ∗k = ⎢
..
⎢
.
⎢
⎢
⎣
O
O

t
L11 L−1∞∞ Σ21
Σ11
Σ21
..
.
O
O

O
Σt21
Σ11
..
.
O
O

···
···
···
..
.
···
···

O
O
O
..
.
Σ11
Σ21

⎤
O
O ⎥
⎥
⎥
O ⎥
⎥
.. ⎥
. ⎥
⎥
⎥
Σt21 ⎦
Σ11となる．これより，第 (1, 2)，(2, 1) 成分において期待される共分散構造と一致しないことがわかる．しかし 3 段階目以降（3 領域目以降）は正確に共分散構造を再現している．つまり X1 , X2 , . . .
の代わりに X3 , X4 , . . . を用いれば良い．
s = 2 の場合 L∗1 の存在が必要になるが，今回は一般的な存在を証明することができなかった（個別の場合に存在を証明する方法は補題 11 を参照）．もしこの存在を仮定すれば
V1 = Σ11 , L11 = Ch(V1 ) ,L∗ki = O
L∗k,k−1 =(1 ≤ i ≤ k − 3) ,lim Lk,k−1 = L∗1 ,k→∞

L∗k,k−2 = Σ31 (Lt∞∞ )−1 = L∗2 ,L∗kk = L∞∞ .

これより領域全体の共分散は(Σ∗k )11 = L11 Lt11 = V1 = Σ11 ,(Σ∗k )21 = L21 Lt11 = Σ21 , (Σ∗k )22 = L21 Lt21 + V2 ,(Σ∗k )31 = L∗2 Lt11 , (Σ∗k )32 = L∗2 Lt21 + L∗1 Lt22 ,(Σ∗k )42 = L∗2 Lt11 ,(Σ∗k )ii = L∗2 (L∗2 )t + L∗1 (L∗1 )t + V∞ = Σ11(i ≥ 3) ,−1(Σ21 − L∗2 (L∗1 )t )t + Σ21 − L∗2 (L∗1 )t = Σ21(Σ∗k )i,i−1 = L∗2 (L∗1 )t + L∗1 Lt∞ = Σ31 V∞(Σ∗k )i,i−2 = Σ31 (Lt∞∞ )−1 Lt∞∞ = Σ31(i ≥ 4) ,(i ≥ 5)となることが分かる．結局 s = 2 の場合，5 領域目以降は期待する共分散構造に従っていることが分かる．つまり X1 , X2 , . . . の代わりに X5 , X6 , . . . を用いれば良い．
補題 11. もし ||Σ31 || ||(Lt∞∞ )−1 ||2 < 1 であれば L∗1 が存在する．
証明．証明は，距離空間 T における写像 Φ : T → T が縮小写像
ある定数 q < 1 があり |Φ(x) − Φ(y)| ≤ q|x − y| ∀x, y ∈ Tブロック行列による正規確率場の逐次シミュレーション法について

143

であれば，
（1）Φ(x∞ ) = x∞ となる不動点 x∞ が一意に存在，（2）任意の点 x0 に対し xn = Φ(xn−1 )によって定義される点列 {xn } は x∞ に収束する，という良く知られた結果に基づく．行列値関数 F (A) = (Σ21 − L∗2 At )(Lt∞∞ )−1 を定義すると||F (A) − F (B)|| = ||L∗2 (At − B t )(Lt∞∞ )−1 ||≤ ||L∗2 || ||A − B|| ||(Lt∞∞ )−1 || = M ||A − B|| ,ここで
M = ||L∗2 || ||(Lt∞∞ )−1 || = ||Σ31 (Lt∞∞ )−1 || ||(Lt∞∞ )−1 || ≤ ||Σ31 || ||(Lt∞∞ )−1 ||2であるから補題が成立する．
6.

数値実験例

今回行ったシミュレーション（表 1，図 2）は，いずれも球型共分散関数を用いた以下の 4 パ
ターンである：
各パラメータの意味は以下の通りである：• dx：隣接点間の水平・垂直距離，• a，b：共分散関数のレンジとシル，• c，r ：各部分領域内の横・縦方向のシミュレーション位置の個数，• n：部分領域数．

（図 3 参照）．V∞ は，行列方程式を既述の繰り返し解法
例 1，2 では a の定義より s = 1 となる
で ||Vn+1 − Vn || ≤ 10−9 となるまで繰り返し求めた．例 3，4 は，共分散関数のレンジがひとつの部分領域の横幅の 2 倍になるため s = 2 となる．例 3，4 では，いずれも Σ31 = O であり，補題 11 から L∗1 は存在する．例えば，例 2 では一つの部分領域中の位置の数は 4, 800 であり，条件付き共分散行列のサイズは 4, 800 × 4, 800 となる．一方で全領域に対する位置の数は 38, 400
で，共分散行列のサイズは 38, 400 × 38, 400 であることに注意する．
次に極限条件付き分散を用いた正当性を確認するため，極限条件付き分散を用いて，各例において発生させたシミュレーション値のヒストグラムと，期待される正規確率密度関数のグラフを重ねたものが，図 4 である．
更に方法 1. 極限条件付き分散を一切用いない，本来の逐次シミュレーション法，方法 2. 行列方程式を解いて極限条件付き分散を予め求め，それを最初から使う，
方法 3. 本来の逐次シミュレーション法で条件付き分散が収束したら，以降はそれを使う，の 3 つの方法での繰り返し回数と計算時間を比較する．方法 1,2 で，極限条件付き分散が求ま
表 1. シミュレーションのパラメータの組み合わせ．

統計数理

144

図 2.

図 3.

第 60 巻

第 1 号 2012

極限条件付き分散を用いたシミュレーション図：上から順に，例 1（21 ∼ 30 領域），例 2（71 ∼ 80 領域），例 3（21 ∼ 30 領域），例 4（71 ∼ 80 領域）．

例 1，2 における領域全体のイメージ図（上），例 3，4 における領域全体のイメージ図（下）．

るまでの繰り返し回数を比較した結果が表 2 である．また，それぞれの例において領域全体をシミュレーションするのにかかった計算時間を，3 つの方法毎に比較した結果が表 3 である．
また c = 80，r = 60 で，領域数と計算時間の関係をプロットしたものが図 5 である．なお今回の数値実験は以下の環境で行った．
（1） CPU：Intel Core i7 CPU 860，
（2）メモり：16 GB，（3） OS：32 bit Lnux，（4）ソフトウェア：R-12.0

ブロック行列による正規確率場の逐次シミュレーション法について

図 4.

各例のシミュレーション値のヒストグラムと平均 0，分散 10 の正規密度関数．例 1（上段左），例 2（上段右），例 3（下段左），例 4（下段右）．

表 2.

条件付き分散の極限を求めるための計算回数の比較．

表 3. 極限条件付き分散の計算方法と計算時間の関係．

145

統計数理

146

図 5.

7.

第 60 巻

第 1 号 2012

例 2 における領域数と計算時間のプロット．横軸は領域数，縦軸はシミュレーション計
算時間（秒）．◦：方法 1，+：方法 2，#：方法 3．

結論

この論文で提案した正規確率場の逐次シミュレーション法は，原理的には任意の領域，点配
置，そして共分散関数に対して正確なシミュレーション法を与える．しかしながら，各段階で登場する行列のサイズを小さくするためには，部分領域とその中の位置集合は適当に小さい必要があり，その分繰り返し回数が多くなる．特に部分領域やその中の位置点集合が合同でな
かったり，共分散関数が有限なレンジを持たぬ場合，毎回条件付き平均・分散を更新する必要がある．
点位置集合が格子状の時に限れば，現在もっとも効率的な方法は離散フーリエ変換と FFT法に基づくアルゴリズムと思われる．しかしこの方法は本質的にトーラス上の確率場に対するアルゴリズムであり，共分散関数が有限でないときは，それを有限範囲の共分散関数で近似する必要があり，近似的な手法になる． FFT 法によるアルゴリズムの計算量は，総位置点の数をn とすると O(n log n) になることが知られている．一方，有限範囲の共分散関数に限れば，今回提案した手法の計算量は図 5 からも推察されるように O(n) になる．但し，実際には毎回かなりの大きさの行列の演算を一定回数必要とするため，この計算量は理論的な興味に留まる．
極限条件付き分散を予め求めておき，それを最初から使う方法は条件付き分散を毎回更新する手間を省き，最初のステップを除けば正確なシミュレーション結果を与えるという点で，今回のアルゴリズムの高速化に有効であると思われるが，対応する行列方程式の解法アルゴリズ
ムは十分な研究が行われていないようで，今回用いた繰り返しアルゴリズムは，各部分領域中の点数が大きくなると，本来の逐次シミュレーションを行い条件付き分散の収束を確認するのと，大差なくなることが分かった．
結論としては，今回提案したアルゴリズムは，離散フーリエ変換法が使える状況では太刀打ちできるものではなく，時間さえかければ実際に実行可能で正確な汎用的逐次シミュレーショ
ン法としての理論的興味という側面が強い．しかしながら条件付きシミュレーション法等では，一度にシミュレーションするには大きすぎる領域上の不規則な配置の点集合上で，多数回

ブロック行列による正規確率場の逐次シミュレーション法について

147

シミュレーションを行う必要がある．そうした場合に使える正確な手法としては，この論文で提案した逐次アルゴリズムが唯一のものと思われる．ここで，今回提案したアルゴリズムは，
各 Xi の成分の幾つかが既知（最初から条件付けられている）の場合でも，適当な変更で有効なことを注意しておく．また，同じ条件で複数回シミュレーションするのであれば，初回で計算した極限条件付き分散や関連した行列を，二回目以降再利用することで高速化が可能になることを最後に注意しておく．

参 考 文 献Benner, P. and Faßbender, H.（2007）. On the solution of the rational matrix equation, EURASIP
Journal on Advances in Signal Processing, 2007, Article ID 21850.
. Geostatistics Modeling Spatial Uncertainty, Wiley & Sons, New
Chilès, J. P. and Delﬁner, P.（1999）
York.
. Production of conditional simulations via the LU triangular decomposition of
Davis, M. W.（1987）
the covariance matrix, Mathematical Geology, 19（2）, 91–98.
. Selected Nonlinear Matrix Equations, VDM Verlag Dr. Müller, Leipzig.
Fital-Akelbek, S.（2009）

ハーヴィル，D. A.（2007）. 『統計のための行列代数上，下』（伊理正夫 監訳）
，シュプリンガー・ジャ
パン，東京．
間瀬茂（2010）. 『地球統計学とクリギング法』，オーム社，東京．
Matheron, G.（1973）. The intrinsic random functions and their applications, Journal of Applied
Probability, 5, 439–468.
. The Fourier integral method: An eﬃcient spectral
Pardo-Iguzquiza, E. and Chica-Olmo, M.（1993）
method for simulation of random ﬁelds, Mathematical Geology, 25（2）, 117–217.
. Matrix Iterative Analysis, Springer Verlag, Heidelberg.
Varga, R. S.（2009）
Wackernagel, H.（2003）. Multivariate Geostatistics, 3rd ed., Springer Verlag, Heidelberg.
Williams, D.（1991）. Probability with Martingales, Cambridge University Press, Cambridge.

148

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

Sequential Simulation of Gaussian Random Fields Based on Block Matrices
Shigeru Mase, Nobuyuki Mochida and Yasuhiro Murayama
Department of Mathematical and Computing Sciences, Tokyo Institute of Technology

As to exact (at least for bounded range covariance function case) simulation methodsof Gaussian random fields, the most direct one is based on Choleski decomposition ofcovariance matrices. But this method becomes intractable as the size of covariance matrix
becomes larger. On the other hand, the most frequently used method based on FFT (Fast
Fourier Transform) can be used only for simulations on regular grids. In this paper, wederive a general sequential formulas for conditional means and covariances of multivariate
Gaussian random vectors and apply them to a sequential (blockwise) simulation methodof Gaussian random fields on general type of locations. The computational cost of thismethod is less than that of FFT method at least theoretically.

Key words: Gaussian random ﬁelds, sequential formulas for conditional means and covariances, block
matrix, sequential simulation.