日本ソフトウェア科学会第 33 回大会 (2016 年度) 講演論文集

統計的モデル検査における偏向分布洗練手法
冨田 尭

素朴な統計的モデル検査では，稀少事象の生起確率を十分な精度で推定する際に，非常に巨大な標本集団を確率的シ
ミュレーションによって生成する必要がある．これを回避するための手法として，偏向確率分布を用いて標本抽出を行う重点抽出法がしばしば用いられており，多くの場合，稀少事象を導く乱数の偏りを基に偏向確率分布を洗練することで効率化が行われている．
本稿では，離散確率分布と指数分布，一様分布，正規分布等の連続確率分布及びその混合によって確率的振る舞いを特徴付けられる一般的な確率的モデルを対象とした偏向確率分布の洗練手法を与える．ここでは特に連続確率分
布に注目し，偏向分布を多変量正規分布で与えることを考える．そして統計的モデル検査のための偏向分布洗練に（変分ベイズ）EM 法を組み込む．また，優良な初期偏向分布を用いなくても洗練手法が機能するように緩和法を導入する．

充足判定ではなくパスに対する充足判定しか必要と1 はじめにしないため，時間有界な線形時間性質であれば一般的

1. 1 背景な（実）時間論理などで記述できないものであっても確率的モデル検査 [4] [13] [10] では，確率的モデルに容易に扱うことができる [17]．さらに，標本抽出と充対して，注目する性質を満たすパスの生起確率を求

足判定は，各標本パスごとに平行に行うことができ，める．統計的モデル検査 [15] [14] [16] [2] [17] は，統計的並列化が非常に容易である．

手法を用いた近似的な有界確率的モデル検査であり，以下の 3 要素からなる．

• 確率的シミュレーションによる標本パス集団の無作為抽出

しかしながら，注目する性質が，それを満たすパスが生起する確率が非常に小さい稀少事象である場合には，その生起確率を十分な精度で推定・検定を行うには非常に巨大な標本集団を生成する必要があり，いか• 各標本パスの充足判定に並列化が容易であったとしても生起確率によっては• 充足確率の統計的推定あるいは仮説検定

現実的な計算資源では対処できないこともあり得る．

数値計算による確率的モデル検査は，無限状態をも重点抽出法（Improtance sampling） [12] はこれつモデルを検査することが困難であり，また，状態爆を回避する有効かつ古典的な手法である．重点抽出法

発といった問題もある．一方で，統計的モデル検査で

では，元の確率分布とは異なり，稀少事象が生起しやは，確率的シミュレーション可能でさえあればどのよ

すい偏向分布を用いて標本を無作為抽出し，確率分布

うなモデルも検査できる [17]．また，モデルに対すると偏向分布の尤度比によって推定結果を補正することで，元の確率分布での稀少事象生起確率を高精度で推Refinement Methods of Biased Distributions for Statistical Model Checking
Takashi Tomita, 北陸先端科学技術大学院大学 情報社会基盤研究センター, Research Center for Advanced Computing Infrastructure, Japan Advanced Institute of
Science and Technology.

定することが可能となる．偏向分布を得るための手法としては交差エントロピー（Cross entropy; CE）
法 [9] [3] [8] がしばしば用いられており，反復により偏向分布を洗練することができる．ただし，ある程度の確率で稀少事象を生起させる優良な初期偏向分布がめの緩和法を導入する．4 章では，提案手法の有効性

与えられないと洗練が正常に進まない．そして，そのを確認するための実験を行い，その結果について考察

ような優良初期偏向分布を得ることは一般には困難

する．そして 5 章で関連研究について述べ，最後に 6である．

でまとめる．

1. 2 内容

2 準備

本研究では，離散確率分布及び連続確率分布に基づ

2. 1 確率過程

いてそのふるまいが規定される確率過程に対して有本稿では，確率過程として，効に機能する偏向分布洗練手法を与える．その離散確• 状態集合 Q = QD ∪ QC ，率分布及び連続確率分布に従う乱数群が確率的シミュ• ふるまいを決定付ける離散確率変数が従う確率

レーションによって得られるパスに一致するため，洗練の果てに得たい偏向分布は，稀少事象が生起しやす
いパスを誘導する乱数群の偏りを説明する分布である．そのような分布を推定することは統計や機械学習[1] の分野で広く行われていることであるが，統計的モ

デル検査での偏向分布洗練では，元の確率分布と偏向
分布の尤度比を考慮する必要がある．離散確率分布部分に関しては既に先行研究 [8] がある．そこで，本研
究では連続確率分布部分に焦点を当てる．確率的モデ

質量関数 D，• ふるまいを決定付ける連続確率変数 x ∈ RK が従う確率密度関数 C，• 現状態の滞在時間関数 S : (QD × Dom(D)) ∪(QC × Dom(C)) → R≥0 ，• 次 状 態 を 与 え る 状 態 遷 移 関 数 T : (QD ×Dom(D)) ∪ (QC × Dom(C)) → Q，
• 初 期 状 態 関 数 I : (Dom(D) → Q) ま た は(Dom(C) → Q)

ルのふるまいを説明する連続確率分布としては，一様からなる 6 つ組 P = ⟨Q, D, C, S, T, I⟩ を扱う．た分布，ガンマ分布（指数分布や Erlang 分布などを含

だし，む），正規分布，Weibull 分布（指数分布や Rayleigh• QD と QC は共通部分をもたない．

分布などを含む）など様々なものがあるが，ここでは• |Dom(D)| は有限であるものとする．

元の連続確率分布がどの種類であっても，その乱数群• C は各次元について独立であるものとする．

の偏りを多変量混合正規分布でとらえることにする．

この確率過程の 1 つの実行は，状態-時間交代列

そして，統計・機械学習で広く用いられている期待q0 t0 q1 t1 · · · で表現でき，これをパスと呼ぶ．有限長

値最大化（Expectation-Maximization; EM）法の状態-時間交代列 ω = q0 t0 q1 t1 · · · qn−1 tn−1 qn を接[5] [1] や変分ベイズ（Variational Bayesian; VB）

頭辞として持つパスの生起確率 Pr P (ω) は以下のよ

EM 法 [1] [6] を，尤度比を考慮して CE 法に組み込ん

うに与えられる．

だ洗練手法を提案する．また，優良な初期偏向分布
を用いずとも洗練が進むように，緩和法を導入する．
そして，提案手法の有効性を確認するため，ごく単純
かつ視覚的にわかりやすい例題を対象とした実験を行った．

1. 3 本稿の構成
2 章では，準備として，本研究で対象とする確率過程の定義を与え，重点抽出法・CE 法の概略を述べる．

3 章では，EM 法・VBEM 法を組み込んだ偏向分布
洗練手法を与え，また，洗練が効率的に進行させるた

Pr P (ω) = Probinit (q0 )
(n−1
)
∏
·
Probstay (qi , ti )
i=0
(n−1
)
∏
trns
·
Prob (qi , qi+1 ) .

(1)

i=0

ただし，

Probinit (q) =

∑


D(d)


d∈{d|I(d)=q}∫
C(x) dx
if I is Dom(D) → Q,
otherwise,{x|I(x)=q}(2)

Probstay (q, t) =

∑


D(d)


d∈{d|S(q,d)=t}∫
C(x) dxさらに，ハイブリッドシステムなど連続的変化をもつ系に関しても，状態 q ∈ Q をロケーションとして

if q ∈ QD ,(3)

その滞在時間中のふるまいの解釈を追加することで対応可能である．

otherwise,{x|S(q,x)=t}
trns
′

Prob

(q, q ) =

∑


D(d)


d∈{d|T (q,d)=q ′ }∫C(x) dx

2. 2 Bernoulli 試行の成功確率推定

if q ∈ QD ,(4)otherwise.

{x|T (q,x)=q ′ }ここで，乱数列 χ = ⟨d0 , x0 ⟩⟨d1 , x1 ⟩ · · · ⟨d2n , x2n ⟩ ∈(Dom(D)Dom(C))n に対して，P 上のパス path(χ) =q0 t0 q1 t1 · · · qn−1
 tn−1 qn は一意に定まる．

I(d0 ) if I is Dom(D) → Q,
q0 =
(5)

I(x0 ) otherwise,S(qi , d2i+1 ) if qi ∈ QD ,
(6)
ti =

S(qi , x2i+1 ) otherwise,T (qi , d2(i+1) ) if qi ∈ QD ,(7)
qi+1 =T (qi , x2(i+1) ) otherwise,
この P の定義では，計算機で確率的シミュレーショ
ン可能な様々な確率過程を表現できる．

• 離散時間 Markov 連鎖– Q が離散，– 任意の q ∈ QD 及び d ∈ Dom(D) について
S(q, d) = 1，
– 任意の q ∈ QC について S(q, x) = 1．
• 連続時間 Markov 連鎖– Q が離散，– 任意の q ∈ QD について T (q, x) ∈ QC かつ S(q, d) = 0，

– 任意の q ∈ QC について T (q, x) ∈ QD か
つ S(q, x) が指数分布に従う（例えば，C のある次元 k が指数分布で S(q, x) = xk ）．

• 離散化された確率微分方程式 [11]（Euler-丸山法）– Q が Euclid 空間（QD が空集合），– C が K 次元標準正規分布，– 任意の q ∈ Q について S(q, x) = ∆t かつ√
T (q, x) = q + F(q) · ∆t + G(Q) · x · ∆t
（F は漂流係数，G は拡散係数）．

Ω を標本空間，f : Ω → R≥0 を Ω 上の確率分布（確率密度関数あるいは関数質量関数），z : Ω → {0, 1}
を Ω 上の指示関数とする．このとき，確率分布 f の下での z(ω) の平均 Ef [z(ω)] 及び分散 Vf [z(ω)] は以下のように与えられる．
∫
Ef [z(ω)] =z(ω) · f (ω) dω,(8)Ω

Vf [z(ω)] = Ef [z(ω)] · (1 − Ef [z(ω)]).

(9)

ω1 , . . . , ωN を確率分布 f に従う独立な N 個の確率変数としたとき，Ef [z(ω)] の Monte Carlo 推定量

γM C は，以下のように与えられる．
N
1 ∑·z(ωn ).
(10)γM C =
N n=1
また，この推定量 γM C の質を示す相対誤差は，以下のように与えられる．

√
V[γM C ]RE [γM C ] =.
(11)
Ef [z(ω)]
ここで，V[γM C ] は以下のように与えられる γM C の分散である．

1· Vf [z(ω)].
(12)
N
そして，独立同分布に従う N 個の ω1 , . . . , ωN を
V[γM C ] =指示関数 z で評価することは，Bernoulli 試行と見なすことができ，その成功確率は平均 Ef [z(ω)] と一致する．ただし，式 (11) より，Ef [z(ω)] ≪ 1 のとき，相対誤差 10−1 程度の成功確率点推定を行うには 100 · (Ef [z(ω)])−1 個の大規模な標本集団が必要となる．
確率的モデル検査では，• 標本空間 Ω はパスの全体集合，• 確率分布 f はパス全体集合上の確率分布，• z(ω) = 1 は ω |= φ（パス ω が線形時間性質 φを満たす）に対応しており，平均 Ef [z(ω)] は φ を満たすパスの∑生起確率を表す．そして，正例数 N
n=1 z(ωn ) は二
項分布に従うため，その生起確率の信頼区間は解析的に求めることができる．

2. 4 交差エントロピー法

2. 3 重点抽出法
前節で述べたように，Ef [z(ω)] ≪ 1 のときに素朴

最適偏向分布に「近い」偏向分布を得るための手法

な統計的モデル検査（平均推定）を行うには，非現実として，交差エントロピー法 [9] がある．交差エント的な規模の標本集団を抽出しなければならない．この

ロピーは以下のように定義され，二つの分布間の「近種の問題を回避する古典的な統計的手法として重点抽出法 [12] がある．z(ω) · f (ω) > 0 ならば g(ω) > 0さ」の評価関数としてしばしば用いられている．
∫
CE (g1 , g2 ) = −g1 (ω) · log g2 (ω) dω(19)である，Ω 上の確率分布 g を考える．このとき，尤

CE 法 で は ，パ ラ メ ー タ 化 さ れ た 偏 向 分 布 の 族

度比関数{gλ | λ ∈ Λ} の中で最適偏向分布 g⋆ との交差エ

f (ω)
.
(13)
L(ω) =g(ω)
に対して，以下が成り立つ．
∫
Eg [z(ω) · L(ω)] =z(ω) · L(ω) · g(ω) dω∫Ω=
z(ω) · f (ω) dω

ントロピー CE (g ⋆ , gλ ) が最も小さい偏向分布 gλ⋆Ω= Ef [z(ω)].

(14)

よって，z(ω) = 1 となるパス ω の生起確率が大きくなるよう偏らせた確率分布（偏向分布）g を用いて重点抽出し，尤度比 L を考慮して補正することで，

Ef [z(ω)] を効率的に推定できる．すなわち，確率分
布 g の下で N 個の標本 ω1 , . . . , ωN を重点的に標本を抽出したとき，母平均 Ef [z(ω)] の重点抽出推定量

γIS は，以下のように与えられる．
N
1 ∑
γIS =
·
L(ωn ) · z(ωn ).
(15)
N n=1
このとき，推定量 γIS の相対誤差 RE [γIS ] は以下のように与えられる．

Ω（を与える最適パラメータ λ⋆ ）を探索する．ここで，最適パラメータ λ⋆ と任意のパラメータ λ′ に関して
以下が成り立つ．

λ⋆ = argmin CE (g⋆ , gλ )λ∈Λ∫z(ω) · f (ω)= argmax· log gλ (ω) dωλ∈Λ ∫Ω Ef [z(ω)]= argmaxz(ω) · f (ω) · log gλ (ω) dωλ∈Λ ∫Ω= argmaxz(ω) · Lλ′ (ω) · gλ′ (ω) · log gλ (ω) dω
λ∈Λ

Ω

(20)
ただし，

f (ω)
(21)
gλ′ (ω)である．従って，あるパラメータ λi に対して，偏向分
Lλ′ (ω) =布 gλi の下で無作為抽出された N 個の標本 ω1 , . . . , ωNがあるとき，最適パラメータ λ⋆ の点推定値 λi+1 は√V[γIS ].
(16)
Ef [z(ω)]
ただし，V[γIS ] は以下のように与えられる，γIS の以下のように与えられる．
N∑z(ωn ) · Lλi (ωn ) · log gλ (ωn )λi+1 = argmax

分散である．
1· (Eg [(z(ω) · L(ω))2 ] − (Ef [z(ω)])2 ).
V[γIS ] =N
(17)
従って，z(ω) · L(ω) の分散が小さい偏向分布 g を用

(22)
このパラメータ更新と更新された偏向分布を用いたRE (γIS ) =いることで効率的な重点抽出が可能となる．なお，重∑み付き正例数 N
n=1 L(ωn ) · z(ωn ) が従う確率分布は一般には未知である．
そして，最適偏向分布は g⋆ は，一般に以下のように与えられる．

z(ω) · f (ω).
(18)
Ef [z(ω)]
ただし，この最適偏向分布は推定しようとしている母g⋆ (ω) =平均 Ef [z(ω)] に依存している．

λ∈Λ

n=1

標本抽出を繰り返すことで，偏向分布を局所最適化できる [3] [8]．
なお，式 (22) は，重み付き対数尤度の最適化問題とみなせる．また，f = gλ となるようなパラメータλ が存在しないようなパラメータ化偏向分布 g を考慮してもよい．

3 提案手法2. 1 節で述べたように，乱数列 χ = ⟨d0 , x0 ⟩⟨d1 , x1 ⟩· · · ⟨d2n , x2n ⟩ ∈ (Dom(D)Dom(C))n に対して，P 上のパス path(χ) = q0 t0 q1 t1 · · · qn−1 tn−1 qn は一意に定まる（式 (5)–(7)）．従って，指示関数 z の台に含まれるパスの生起確率 PrP ({ω | z(ω) = 1}) は，そのxn,1 , . . . , xn,Mn とする．

ようなパスを導出する乱数列の生起確率と一致する．

元の連続確率分布 C の定義域 Dom(C) が K-次元

この乱数列の生起確率は，離散確率質量の総積と連続

Euclid 空間全体 RK でない場合でも，通常の多変量

確率密度の総積の積の形をもっており，式 (22) より，

正規分布を用いる．ただし，無作為抽出や尤度比（の離散分布 D と連続分布 C がそれぞれ互いに独立なパ中の確率密度）計算の際には，平均と分散共分散行列

ラメータセット λD と λC でパラメータ化することは同じで定義域を Dom(C) に制限した切断多変量正
規分布 [7] として扱う．ある程度最適偏向分布に近けで，それぞれ平行に最適化可能である．
離散分布 D のパラメータ化及び CE 法を用いた偏向分布洗練手法は Jegourel らの手法 [8] などを用い

れば偏向分布でも重点抽出法に十分利用できるため，
この近似を用いる．

ればよいため，以降では，連続分布 C のパラメータ

なお，C の次元数が K であることから，M は基本化及び偏向分布洗練手法を考える．既に述べたよう的には K 変量混合正規分布である．ただし，実用上，に，式 (22) より，乱数列 χ に現れる連続分布乱数

重み付き対数尤度に影響を与えないあるいは影響を

x0 , x1 , . . . , x2n をいくつかの群に分けてそれぞれ独与えるにしてもそれが非常に小さい次元に関しては

立にパラメータ化し，平行に偏向分布の最適化を行う

捨象してもよい．

ことができる．乱数群の分け方としては，• 滞在時間の決定に用いられたかどうか，

3. 1 期待値最大化法を利用した偏向分布洗練• 次状態の決定に用いられたかどうか，

観測データに対する対数尤度の最適化問題を解く基• 滞在時間/次状態の決定に用いた際の現状態がど本的な手法として，期待値最大化法（EM 法） [5] [1]の状態群（⊆ QC ）に含まれているか，がある．重点抽出法を用いた統計的モデル検査のため

あるいは，それらの組み合わせ等が考えられる．

の多変量混合正規偏向分布洗練に EM 法を組み込ん

本稿では，各乱数群の偏りを捉えるための分布として，多変量混合正規分布 M を用いることとする．
L∑
Mλ (x) =πl · Nµl ,Σl (x),(23)

手続き 1. 与えられた初期偏向分布パラメータ λinit
を現在の偏向分布パラメータ λold として初期化し，パラメータ収束条件を満たすまで次の (1)–(3) を繰り

l=1

λ = {πl , µl , Σl }1≤i≤L .

だ手法として，以下の手続きを与える．

(24)

こ こ で ，L は 混 合 要 素 数 ，π1 , . . . , πL は 混 合 率∑（ 1≤l≤L πl = 1），Nµ,Σ は平均が µ かつ分散共分散行列が Σ の多変量正規分布である．

Nµ,Σ (x) =()(x − µ)T · Σ−1 · (x − µ)1√.
· exp −
2(2π)n · det(Σ)(25)
つ ま り，偏 向 分 布 の 洗 練 は ，抽 出 さ れ る 標 本

返す．

1. 現在の偏向分布パラメータ λold = {πlold , µoldl ,
Σoldl }1≤l≤L で与えられる偏向分布 Mλold を用いて N 個の標本パス ω1 , . . . , ωN を重点抽出する．

2. 各標本パス ωn 中の注目する各乱数標本 xn,1 ,. . . , xn,Mn に対して混合要素の負担率 rl,n,m を

パ ス ω1 , . . . , ωN に 対 す る 重 み 付 き 対 数 尤 度求める（E ステップ）．
πlold · Nµold ,Σold (xn,m )llrl,n,m =(27)
Mλold (xn,m )3. E ステップで求めた負担率から，偏向分布パラ

θ(λ, {ω1 , . . . , ωN }) を最大化する多変量混合正規分布

メータを λnew = {πlnew , µnew
, Σnew}1≤l≤L に更ll（のパラメータ λ）を求めることに相当する．

新する（M ステップ）．
N∑
Lλold (ωn ) · Rl,nθ(λ, {ω1 , . . . , ωN }) =∑z(ωn ) · Lλ (ωn ) ·
1≤n≤N

∑

log(Mλ (xn,m )).

1≤m≤Mn(26)
ここで，標本パス ωn 上に現れる注目する乱数群を

πlnew =

n=1
N
∑
n=1

,
Lλold (ωn ) · Mn

(28)

N
∑

µnew
=
l

n=1

Mn
∑

3. 2 変分ベイズ EM 法を利用した偏向分布洗練

rl,n,m · xn,m

m=1
N
∑

N
∑

=
Σnew
l

Lλold (ωn ) ·対数尤度の最適化問題を解くもう一つの基本的な手

,

(29)

Lλold (ωn ) · Rl,nn=1

Lλold (ωn ) ·

n=1

Mn
∑

ある．多変量正規分布 Mλ においては，混合率 π の共役事前分布は Dirichlet 分布，各混合要素の平均⊗
rl,n,m · (xm − µold
l )

m=1
N
∑

.

Lλold (ωn ) · Rl,n
(30)
Mn
∑

µl と精度行列 Σ−1の共役事前分布は正規-Wishartl分布である．

EM 法と同様に，重点抽出法を用いた統計的モデル

n=1

ただし，

法として，変分ベイズ EM 法（VBEM 法） [1] [6] が検査のための，多変量混合正規偏向分布洗練に VBEM
法を組み込んだ手法として，以下の手続きを与える．

(31)

手続き 2. 初期事前分布パラメータを {αlinit , βlinit ,Lλold (ωn ) = z(ωn ) · Lλold (ωn ),(32)

νlinit , ξlinit , Wlinit }1≤l≤L とする．与えられた初期偏向y⊗ = y · (yT ).

(33)

分布パラメータ λinit を現在の偏向分布パラメータ■λold 初期事前分布パラメータの {ξlinit , Wlinit }1≤l≤LRl,n =

rl,n,m ,

m=1

式 (22) より，標本 xn,m が Lλold (ωn ) 個観測されを現在の事前分布パラメータ {ξlold , Wlold }1≤l≤L とたものとして解釈できるため，式 (28)–(30) は重み付して初期化し，パラメータ収束条件を満たすまで次のきの更新式となっている．また，式 (30) で µold
を用l(1)–(3) を繰り返す．

いているため，各標本パス ωn に対して，式 (27) 及

1. 現在の偏向分布パラメータ λold = {πlold , µoldl ,

び式 (28)–(30) の分子・分母の途中計算を独立に行う

Σoldl }1≤l≤L で与えられる偏向分布 Mλold を用いことができる．そのため，通常の EM 法と異なり，標て N 個の標本パス ω1 , . . . , ωN を重点抽出する．

本集団全体の情報を保持する必要がなく，標本抽出及

2. 各標本パス ωn 中の注目する各乱数標本 xn,1 ,び E ステップの計算，M ステップの途中計算を並列′. . . , xn,Mn に対して混合要素の負担率 rl,n,mを化できる．そして，通常の EM 法では固定された観測データに対して EM サイクルを回すが，反復が進むごとにより洗練された標本パス集団を用いる方が効率的なため，この手続きでは 1EM サイクルを 1CE
サイクルとしている．
通常の EM 法と同様に，式 (22) で与えられる重み付き対数尤度 θ の変化を観察することで，パラメータの収束を確認できる．ただし，反復ごとに異なる標本パス集団を用いて洗練を行うため，対数尤度は単調増加しない．また，パラメータの更新に µoldを用いlているため，収束には通常より多くの反復数が見込ま求める（E ステップ）．
exp(ρl,n,m )′.
rl,n,m= L∑
exp(ρl′ ,n,m )
l′ =1

ただし，

( L
)
∑ initρl,n,m = ψ(αlinit ) − ψαl′( K (l′ =1)∑1
νlinit + 1 − k
+ ·
ψ
2
2
k=1

+ log(det(Wlold ))+ K · ((βlinit )−1 − ln(π))

味がある標本パスは，z(ω) = 1 となる正例 ω のみである．そのため，ステップ 1 の重点抽出においてはN 個の正例が得られるまで抽出を繰り返すことでより効率化できる．

)− νlinit · M(xn,m , ξlold , Wlold )れる．
そして，式 (22) より，パラメータ更新において意

(34)

.
(35)

ここで，ψ はディガンマ関数である．
(
)
u
1 ∑ 1
ψ(s) = lim log(u) − −.
u→∞s v=1 s + v
また，M は以下のように定義される．

(36)

M(x, y, W) = (x − y)T · W · (x − y).

(37)

3. E ステップで求めた負担率から，偏向分布パラ

メータを λnew = {πlnew , µnew
, Σnew}1≤l≤L に，ll

標本数に比例する．そのため，EM 法の場合と同様に

事前分布パラメータを {ξlnew , Wlnew }1≤l≤L に更

1EM サイクルが 1CE サイクルとなっており反復が新する（M ステップ）．
N∑′
Lλold (ωn ) · Rl,n

進むごとにより洗練に適した標本パス集団を用いる

n=1
πlnew = N∑,

Lλold (ωn ) ·

n=1
µnew
=
l

N
∑

Σnew
=
l

(38)

Mn
∑

次の反復で利用できるようにスケーリングしている．
′
rl,n,m· xn,mm=1
N∑

期化している（つまり更新していない）．また，それに合わせて，項 (νlinit + Nl′ )/νlinit によって Wlnew を

Lλold (ωn ) · Mn

n=1

N
∑

この手続きでは，これらのパラメータは反復ごとに初この手続きでは，式 (38)–(40) で与えられる観測

,

(39)

データに対する点推定値によって偏向分布パラメー′
Lλold (ωn ) · Rl,nタ λ を更新している．N → ∞ のとき，この点推定

Mn∑

値は偏向分布パラメータの期待値に収束する．また，

n=1

Lλold (ωn ) ·

n=1

′
⊗
rl,n,m· (xm − µold
l )

m=1
N
∑

パラメータの収束は，変分下界や式 (26) の重み付き,′
Lλold (ωn ) · Rl,n

対数尤度の変化を観察することで判断できる．ただし，EM 法の場合と同様に，これらはやはり単調増加

n=1

(40)

しない．そして，反復ごとに異なる標本パス集団を用β init · ξlold + Nl′ · µnewlξlnew = l,(41)
βlinit + Nl′
init
′ (
ν
+N
Wlnew = l init l · (Wlold )−1 + Nl′ · Σnewlνl)−1
βlinit · Nl′+ init· (µnew− ξlold )⊗lβl + Nl′

いていることから，事前分布パラメータは適宜初期化(42)

た洗練手法を与えた．しかしながら，初期偏向分布が

ただし，
′
Rl,n
=

Mn
∑
m=1
N
∑

Nl′ =

(43)

·
Lλold (ωn )3. 1–3. 2 節では，与えられた初期偏向分布を元にし手続きは有効に機能しない．そのため，優良な初期偏向分布を零から構成する手法，あるいは，容易に構成′
Lλold (ωn ) · Rl,nN∑

3. 3 緩和法

正例を生起させる確率が非常に小さい場合，これらの′rl,n,m,n=1

してもよい．

N∑z(ωn ).

(44)

仕組みが必要である．

n=1

n=1

可能な分布（あるいは元の分布）からでも洗練が進む■

なお，α· は事前 Dirichlet 分布の集中度パラメータであり，l 番目の混合要素に関して，νl· （resp. Wl· ）は事前 Wishart 分布の自由度（resp. 尺度行列），βl·
は尺度パラメータ，ξl· は事前正規分布の平均である．
基本的な特徴は，手続き 1 と同様である．

Nl′ は負担正例数を示す．通常の VBEM 法とは異なり重み Lλold (ωn ) を考慮する必要があり，稀少事象を捉えるパラメータ λ に対してはこの重みは非常に小さい．そのため，事前観察標本数とも解釈できる αなどがこの重みと比べて小さくならないようスケーリ
ングを行っている．また，事前分布パラメータの内，α，β ，ν の更新値は，本来，各混合要素ごとの負担

そこで本稿では，容易に構成可能な分布からでも洗練が進むように，緩和法を導入する．基本的なアイ
ディアは• 洗練の途中段階では，対数尤度が高い偏向分布
ではなく，z の台に近いパスが生起しやすい偏向分布が得られやすいように重み付けを行うことである．具体的なアプローチは以下の通りである．

• 指示関数 z （= z0′ とする）の代わりに以下のよ′うな特徴をもつ評価関数 z1′ , . . . , zH
: Ω → [0, 1]を用いる．

– 任意の h ∈ {1, . . . , H} について，′
Supp(zh−1) ⊂ Supp(zh′ ).

(45)

ただし，Supp は台を与える関数である．

たは 2 のステップ 1 の標本抽出時に一定の割合以上

Supp(zh′ ) = {x ∈ Dom(zh′ ) | f (zh′ ) > 0}.
(46)′– 任意の h ∈ {1, . . . , H}，x ∈ Supp(zh−1)，′
y ̸∈ Supp(zh−1) について，′の標本パスが zh−1の台に乗った次の反復のときなどが考えられる．
なお，C が多変量正規分布でない場合，混合数として適当な L を定め，C に従う乱数を適当数無作為抽

zh′ (x) ≥ zh′ (y).

(47)′
– 元の分布 f において，Supp(zH) に含まれるパスの生起確率は小さくない．

• まず，初期偏向分布を元の分布 C，h := H として，洗練の進み具合に応じて h を小さくしなが出し，それを k 平均法などで L 個のクラスタに分け，各クラスタのサイズから混合率を，各クラスタの標本
平均と標本分散共分散から各混合要素の平均と分散
共分散行列を決め初期偏向分布パラメータとすればよい．

ら（緩和を弱めながら）評価関数 zh′ を z の代わりに用いて偏向分布洗練（手続き 1 または 2）を
行う．

4 実験
3 節の提案手法を実装し，有効性を示すための実験• 洗練の途中段階では，z の台に乗る（あるいはそを行った．

れにより近い）パスが生起しやすい偏向分布を得ることを優先し，式 (32) 中の尤度比 Lλold (ωn )4. 1 実装及び実験環境を 1 で置き換える．

提案手法の実装には R を用いた．多変量正規分布• 最終的に，十分高い確率で z （= z0′ ）の台に乗の確率密度計算と無作為抽出には mvtnorm ライブラるパスが生起するようになった時に得られている

リ，ディガンマ関数の計算や各種行列演算には R 標

偏向分布を初期分布として緩和なしの偏向分布

準の関数を利用した．そして，元の分布から無作為洗練を行う．

抽出した標本集団を R 標準の k 平均クラスタリング

緩和評価関数 zh′ の設定の仕方には様々のものが考関数 kmeans を利用して分類し，クラスタのサイズ

えられる．例えば，元の指示関数 z が線形時間論理

比を初期偏向分布の混合率 π init ，各クラスタの標本式 φ で特徴づけられるものである場合，その過大近平均及び標本分散共分散を各混合要素の平均 µinitとl

似系列 φ1 , . . . , φH を用いて以下のように与えること

分散共分散 Σinitとして初期化している．各事前分布
l

ができる．


1




′
zh (ω) = 0min{h′ | ω |= φh′ }


′ ≤h
1≤h

1 −
h+1

パラメータの規定値は，各混合要素に対して αlinit は

if ω |= φ,if ω ̸|= φh ,otherwise.

(48)
ただし，h が小さい φh ほど φ に近い性質とする．
途中段階でも尤度比 Lλold を考慮する場合，最適偏向分布に近いほど尤度比が小さくなるため，緩和評価
関数 zh′ の値が大きい標本でも尤度比全体としての値が大きくなるとは限らない．そのような場合，z の台に乗る標本が高確率で得られる分布まで洗練が進まない可能性がある．
評価関数 zh′ を用いて偏向分布洗練を行っている際の緩和を弱めるタイミングに関しては，手続き 1 ま

10−3 ，βlinit は 10−1 ，νlinit は K ，ξlinit は零ベクトル，
Wlinit は単位行列とした．洗練の途中段階では 1EM
サイクルで利用する標本パスは正例 100 個とした．また，洗練の最終段階（緩和なし）では，1EM サイク
ルで利用する標本パスは，正例 10000 個または負例

10000 個が観測されるまでに得られた正例とし，1 個も正例が観測できなかった場合には 1 サイクル前のパラメータセットに戻している．そして，同一の評価関数を用いている間を 1 段階とし，各段階に関して300EM サイクルを上限とした．現在の評価関数が zh′（h > 0）のとき，標本抽出時に 95%以上の標本パス′が zh−1の台に乗ったときにはその際に利用したパラ

メータセットを途中解として次の段階に移行する（評′価関数を zh−1に変更する）．また，標本パス集団に対する（標本サイズあたりの）対数尤度の最大値が 5サイクル以上更新されなかった場合，パラメータは十分収束したものとみなし，最大値を与えた時のパラ
メータセットを最適パラメータとして次の段階に移行する（あるいは解として出力する）．なお，段階を進める際には，事前分布パラメータ ξl と Wl を初期化するようにした．
実験は，OS X El Capitan (10.11.6)，2.6GHz Intel

Core i5 CPU，8GB 1600MHz DDR3 メモリの Mac
mini (Late2014) 計算機上で行った．
4. 2 実験内容
提案手法を利用することで優良な初期偏向分布がなくても偏向分布洗練が正常に進むことを確認する
ため，ごく単純で視覚的にわかりやすい例を用いる．

図1

分布 Mλorg に従う無作為抽出標本（10000 個）の確率過程のふるまいを規定する連続確率分布 C は

散布図

以下のパラメータ λorg = {π org , µorg , Σorg } をもつ
混合数 4 の 2 変量正規分布 Mλorg とする．

= ⟨0.1, 0.2, 0.3, 0.4⟩,
⟨[
]T [
]T [
µorg =
10 5 , 7 −1 , 5

以上のような指示関数 z（及び緩和評価関数 zh′ ）と

org

(49)
連続確率分布 C = Mλorg に対して，EM 法に基づい]T ⟩, た手続き 1 と VBEM 法に基づいた手続き 2 のそれぞ0
れについて，C を初期偏向分布として，
(50)
 
 
 
⟩
⟨
• 途中段階では尤度比 Lλold を無視し，重みを zh′
2 0
3 2
1/2
0
5 1
org
,
,
,
 .
Σ = 
とする強い緩和を適用した場合（実験 (i)）と
0 2
2 3
0
1/2
1 5
• 途中段階でも尤度比 Lλold を考慮し，重みを(51)zh′ · Lλorg とする弱い緩和を適用した場合（実験
この混合正規分布に従う標本の無作為抽出散布図を(ii)）
図 1 に示す．
π

]T [
6 , 0

また，議論を単純化するため，パス ω は Mλorg にについて偏向分布洗練を行い，その偏向分布を用いて従う 10 個の乱数 x1 , . . . , x10 のみによって特徴づけ

標本数 100000 の重点抽出平均推定を行う実験を実施られるものとする．そして，指示関数 z（= z0′ ）及びした．

′
緩和評価関数 zh を以下のように設定する．
1
if d(ω) ≤ 10,′zh (ω) = 0
if d(ω) > 10 · (h + 1),
d(ω)/10 − 1otherwise.
1−
h
(52)
∑10
′
ただし，d(ω) = i=1 |xi − [6 3]|. 直観的には，zh (ω)

は，10 個の確率変数の点 (6, 3) までの平均距離が 1以下ならば 1 であり，平均距離が 1 から離れるほど線形に減少し，1 + h 以上離れた時には 0 となる．なお，点 (6, 3) は分布 C の谷の部分にあたる．

4. 3 実験結果
実験 (i) 及び実験 (ii) において，各段階の計算に要した時間（秒）と反復数をそれぞれ表 1 と表 2 に示す．なお，反復数にアスタリスク ∗ が付与されている
項目は，パラメータの収束ではなく，95%以上の標本パスが次段階の評価関数の台に乗ったことによる段階進行を表す．また，参考として，各段階で得られた偏向分布に従う標本の無作為抽出散布図を付録 A に示す．
表 1 と表 2 より，EM 法に基づいた手続き 1 でも

表1

途中段階で尤度比を無視する強い緩和を用いた場合（実験 (i)）の各段階の計算に要した時間及び反復数
重み

EM
実行時間
反復数

実験 (i) 及び実験 (ii) の結果として得られた偏向分布を用いた重点抽出平均推定（標本数 100000）の結果

VBEM
実行時間
反復数

z6′
z5′
z4′
z3′
z2′
z1′

3.920
2.789
6.660
7.758
8.804
8.836

4*
3*
7*
8*
9
9

8.688
10.47
10.32
15.64
15.45
17.13

5*
6*
6*
9
9
10

z · Lλold

2162

13

1103

9

表2

表3

実行時間

推定平均

相対誤差

正例率

EM1
EM0

885.3
886.5

2.98e-22
2.74e-22

.0544
.0137

.0287
.5496

EM’1
EM’0

897.1
900.3

2.13e-22
2.73e-22

.3634
.0129

.0001
.5419

VBEM1
VBEM0

893.3
891.6

2.55e-22
2.73e-22

.0523
.0314

.1789
.4325

VBEM’1
VBEM’0

887.5
905.9

2.12e-22
2.89e-22

.3865
.0337

.0004
.4443

途中段階でも尤度比を無視しない弱い緩和を用いた

場合（実験 (ii)）の各段階の計算に要した時間及び反復数

いても，途中段階で尤度比を無視しても考慮しても最終的な偏向分布の質は同程度であるといえる．尤

重み

EM’
実行時間
反復数

VBEM’
実行時間
反復数

z6′ · Lλold
z5′ · Lλold
z4′ · Lλold
z3′ · Lλold
z2′ · Lλold
z1′ · Lλold

10.65
15.17
15.34
36.44
18.00
78.09

6*
9
9
21
9
15

21.69
17.76
22.74
23.11
19.28
22.60

12
9
12
12
9
12

z · Lλold

2308

14

1675

14

度比を無視した場合と考慮した場合の大きな違いは，最終段階に入った際の偏向分布の正例率である．今回の実験では洗練が正常に進んだが，尤度比を考慮した場合には段階を進めた（評価関数の台が狭まった）直後の反復で正例が抽出できずに洗練が停止する危険性を示唆している．以上のことから，途中段階で尤度比を無視することは，洗練をより確実に進行させるという点で有効であると考えられる．

VBEM 法に基づいた手続き 2 でも，尤度比を無視す
る強い緩和を用いることで途中段階の洗練に要する5 関連研究時間及び反復数が多少削減されていることがわかる．

交差エントロピー法の統計的モデル検査への適用

ただし，計算コストの大部分は最終段階の緩和なしのは Clarke と Zuliani によって提案された [3]．彼らは，部分であるため，計算時間という面では尤度比を考慮

解析的に解くことのできる指数分布族の確率分布をしないというヒューリスティックの有効性は低い．な対象とした手法を与えている．一方，本稿では，指数

お，尤度比を考慮する場合において反復数に対して実分布族には含まれない多変量正規分布を偏向分布と行時間が長いのは，負例の生起確率が高く，100 個のして扱っている．

正例を得るまでに多くの標本パスを必要としているためである．

離散時間 Markov 連鎖に対する統計的モデル検査への適用は Jegourel らによって提案されている [8]．

そして，EM 法に基づいた手続き 1 及び VBEM

離散時間 Markov 連鎖のふるまいは離散確率分布の法に基づいた手続き 2 の，途中段階で尤度比を無視

みによって規定されるものであり，彼らの手法は離した場合（EM と VBEM）と考慮した場合（EM’ と散確率分布を対象としている．それに対し，本稿で

VBEM’）それぞれについて，緩和なしの最終段階には，任意の連続確率分布を対象とした手法を与えた．

入った際の偏向分布を用いた重点抽出平均推定（EM1，

Jegourel らの手法と組み合わせることで，離散確率

VBEM1，EM’1，VBEM’1）と最終的に得られた偏分布と連続確率分布のどちらによっても振る舞いが規

向分布を用いた重点抽出平均推定（EM0，VBEM0，

定される確率過程も扱うことができる．

EM’0，VBEM’0）の結果を表 3 に示す．
表 3 より，EM 法と VBEM 法のどちらの場合につまた，本稿では，先行研究 [3] [8] では与えられていなかった初期偏向分布問題への対策として緩和法も提

案した．

6 まとめ
本稿では，稀少事象の生起確率推定に用いられる重点抽出法に利用する偏向分布を洗練手法を提案した．
離散確率分布と連続確率分布によってそのふるまいが規定される様々な確率過程を扱うことができる．連続確率分布部分に対しては，偏向分布として多変量正規分布を用い，EM 法・VBEM 法という基本的な手法を，統計的モデル検査の並列化容易性を損なわないように改変して反復 CE 法に組み込んだ．また，優良な初期偏向分布を与えなくても洗練が進むように緩和法を導入した．これにより，洗練の途中段階では，稀少事象か否かを与える指示関数を稀少事象に近いほど高い値を与える評価関数に変更し，また，尤度比を考慮しないことで，稀少事象が生起しやすい偏向
分布を効率的に得ることができる．そして，実験によりその有効性を確認した．

参 考 文 献[ 1 ] Bishop, C. M.: Pattern Recognition and Machine Learning, Springer-Verlag New York, 2006.
[ 2 ] Bulychev, P., David, A., Larsen, K. G., Legay,
A., Li, G., Bøgsted Poulsen, D., and Stainer,
A.: Monitor-Based Statistical Model Checking for
Weighted Metric Temporal Logic, Logic for Programming, Artificial Intelligence, and Reasoning,
Bjørner, N. and Voronkov, A.(eds.), Lecture Notes
in Computer Science, Vol. 7180, Springer Berlin
Heidelberg, 2012, pp. 168–182.
[ 3 ] Clarke, E. M. and Zuliani, P.: Statistical Model
Checking for Cyber-Physical Systems, AutomatedTechnology for Verification and Analysis, Bultan,
T. and Hsiung, P.-A.(eds.), Lecture Notes in Computer Science, Vol. 6996, Springer Berlin Heidelberg, 2011, pp. 1–12.
[ 4 ] Courcoubetis, C. and Yannakakis, M.: The
Complexity of Probabilistic Verification, J. ACM,
Vol. 42, No. 4(1995), pp. 857–907.
[ 5 ] Dempster, A. P., Laird, N. M., and Rubin,
D. B.: Maximum Likelihood from Incomplete Data
via the EM Algorithm, Journal of the Royal Statistical Society. Series B (Methodological), Vol. 39,
No. 1(1977), pp. 1–38.
[ 6 ] Fox, C. W. and Roberts, S. J.: A tutorial on variational Bayesian inference, Artificial Intelligence
Review, Vol. 38, No. 2(2012), pp. 85–95.
[ 7 ] Genz, A. and Bretz, F.: Computation of Multi-

variate Normal and t Probabilities, Lecture Notes
in Statistics, Vol. 195, Springer-Verlag Berlin Heidelberg, 2009.
[ 8 ] Jegourel, C., Legay, A., and Sedwards, S.:
Cross-Entropy Optimisation of Importance Sampling Parameters for Statistical Model Checking,Computer Aided Verification, Madhusudan, P. andSeshia, S.(eds.), Lecture Notes in Computer Science, Vol. 7358, Springer Berlin Heidelberg, 2012,
pp. 327–342.
[ 9 ] Kroese, D. P., Rubinstein, R. Y., Cohen, I.,
Porotsky, S., and Taimre, T.: Cross-Entropy
Method, Springer US, 2013, pp. 326–333.
[10] Kwiatkowska, M., Norman, G., and Parker, D.:
Stochastic Model Checking, Formal Methods for
Performance Evaluation, Bernardo, M. and Hillston, J.(eds.), Lecture Notes in Computer Science, Vol. 4486, Springer Berlin Heidelberg, 2007,
pp. 220–270.
[11] Øksendal, B.: Stochastic Diﬀerential Equations:
An Introduction with Applications, Springer-Verlag
Berlin Heidelberg, 2003.
[12] Srinivasan, R.: Importance Sampling: Applications in Communications and Detection, SpringerVerlag Berlin Heidelberg, 2002.
[13] Vardi, M.: Probabilistic Linear-Time Model
Checking: An Overview of the Automata-Theoretic
Approach, Formal Methods for Real-Time andProbabilistic Systems, Katoen, J.-P.(ed.), Lecture
Notes in Computer Science, Vol. 1601, Springer
Berlin Heidelberg, 1999, pp. 265–276.
[14] Younes, H. L., Kwiatkowska, M., Norman,
G., and Parker, D.: Numerical vs. Statistical
Probabilistic Model Checking, International Journal on Software Tools for Technology Transfer,Vol. 8(2006), pp. 216–228.
[15] Younes, H. L. and Simmons, R. G.: Probabilistic Verification of Discrete Event Systems Using Acceptance Sampling, Computer Aided Verification, Brinksma, E. and Larsen, K. G.(eds.), Lecture Notes in Computer Science, Vol. 2404, Springer
Berlin Heidelberg, 2002, pp. 223–235.
[16] Zuliani, P., Platzer, A., and Clarke, E. M.:
Bayesian Statistical Model Checking with Application to Simulink/Stateflow Verification, Proceedingsof the 13th ACM International Conference on Hybrid Systems: Computation and Control, HSCC ’10,
New York, NY, USA, ACM, 2010, pp. 243–252.
[17] 冨田尭, 萩原茂樹, 伊藤宗平, 米崎直樹: 確率頻度時
間論理の統計的モデル検査, コンピュータソフトウェア,
Vol. 31, No. 3(2014), pp. 336–356.

A 実験の各段階で得られた偏向分布
参考として，4 章の実験 (i) 及び実験 (ii) の各段階で得られた偏向分布に従う標本の無作為抽出散布図を図 2-5 示す．

(a) 第 1 段階（z6′ ）(b) 第 2 段階（z5′ ）(e) 第 5 段階（z2′ ）図2(c) 第 3 段階（z4′ ）(f) 第 6 段階（z1′ ）(d) 第 4 段階（z3′ ）(g) 最終段階（z · Lλold ）EM 法に基づいた手続き 1（途中段階では尤度比を無視）の各段階で得られた偏向分布に従う無作為抽出標本（1000 個）の散布図(a) 第 1 段階（z6′ · Lλold ） (b) 第 2 段階（z5′ · Lλold ） (c) 第 3 段階（z4′ · Lλold ） (d) 第 4 段階（z3′ · Lλold ）(e) 第 5 段階（z2′ · Lλold ） (f) 第 6 段階（z1′ · Lλold ）図3(g) 最終段階（z · Lλold ）EM 法に基づいた手続き 1（途中段階でも尤度比を考慮）の各段階で得られた偏向分布に従う無作為抽出標本（1000 個）の散布図

(a) 第 1 段階（z6′ ）(b) 第 2 段階（z5′ ）(e) 第 5 段階（z2′ ）図4(c) 第 3 段階（z4′ ）(f) 第 6 段階（z1′ ）(d) 第 4 段階（z3′ ）(g) 最終段階（z · Lλold ）VBEM 法に基づいた手続き 1（途中段階では尤度比を無視）の各段階で得られた偏向分布に従う無作為抽出標本（1000 個）の散布図(a) 第 1 段階（z6′ · Lλold ） (b) 第 2 段階（z5′ · Lλold ） (c) 第 3 段階（z4′ · Lλold ） (d) 第 4 段階（z3′ · Lλold ）(e) 第 5 段階（z2′ · Lλold ） (f) 第 6 段階（z1′ · Lλold ）図5(g) 最終段階（z · Lλold ）VBEM 法に基づいた手続き 1（途中段階でも尤度比を考慮）の各段階で得られた偏向分布に従う無作為抽出標本（1000 個）の散布図