PRISMは生まれたての計算機言語であるので、現段階ではPRISMに関する外部発 表は行なっていない。
本研究で試作したPRISM処理系を用いて前述した血液型遺伝子モデル、隠れマ
ルコフモデル、ベイジアンネットワーク、回路故障モデルという4つの統計モ
デルを対象にモデルの記述、確率計算、確率パラメータ学習について実験を行
なった。その結果HMMからベイジアンネットワークに至る広範囲の統計モデル
がPRISMで記述でき、PRISM処理系(学習系)により確率パラメータ学習も可能で
あることが確認された。以下に各々の統計モデルに関する実験結果をまとめる。
(1) 血液型遺伝子モデル
(2)BSプログラムの記述例 の血液型遺伝子モデルに制御宣言とユーティリティ
プログラムを加えた PRISM プログラムを以下に示す。ターゲット宣言(
target/2)によりターゲット述語は bloodtype/1 と指定され、データ
宣言( data/1)により教師データファイルは japanese.datと指定
される。
% Blood type model: bloodtype(a) :- genotype(a,a); genotype(a,o); genotype(o,a). bloodtype(b) :- genotype(b,b); genotype(b,o); genotype(o,b). bloodtype(o) :- genotype(o,o). bloodtype(ab) :- genotype(a,b), genotype(b,a). genotype(X,Y) :- gene(father,X), gene(mother,Y). gene(P,a) :- bsw(1,P,1). gene(P,b) :- bsw(1,P,0), bsw(2,P,1). gene(P,o) :- bsw(1,P,0), bsw(2,P,0). % Control declarations: target(bloodtype,1). data('japanese.dat'). % Utility program: print :- prob(bloodtype(a),Pa), write('P(A) = '), write(Pa),nl, prob(bloodtype(b),Pb), write('P(B) = '), write(Pb),nl, prob(bloodtype(o),Po), write('P(O) = '), write(Po),nl, prob(bloodtype(ab),Pab), write('P(AB)= '), write(Pab),nl. sampler(0,[]). sampler(N,[X|Samples]) :- N>0, sample(bloodtype(X)), N1 is N-1, sampler(N1,Samples).
そして教師データファイル japanese.dat に bloodtype(・) を日本人の血液型の比率と同じく bloodtype(a):bloodtype(b):bloodtype(o):bloodtype(ab)=4:2:3:1 の割合で格納し、 gene(p,・) の確率分布を推定す る。
学習系で推定を行なうとパラメータ値はおよそ5回の更新で収束した。学習結 果としてパラメータ(グループ ID は 1, 2)の収束値 Conv とそのパラメータの下での確率分布 PDB(gene(p,・)=1),PDB(bloodtype(・)=1) を表 --2 に、学習中のゴールの連言の対数尤度(logarithmic likelihood)の推移グラフを 図 --5 に示す。横軸が EMの繰り返しステップ数、縦軸が対数尤度である。この学習結果より日本人全 体では血液型遺伝子 A, B, O の分布はおよそ 0.292:0.163:0.545 であると 推定される。推定された遺伝子の分布から表現型(A, B, O, AB)の分布を計算 すると教師データに近い比が得られる。
またこのプログラムではユーティリティプログラムとして print/1 と sampler/2 を定義している。 print/1 は血液型の確率分布を標準 出力に出力するものである。また sampler/2 は sampler( +N, -Samples) の形で与えられたとき、血液型を N 個サン プリングしたものを Samples に返す。例えば次のように質問できる。
| ?- print. P(A) = 0.40387451889260306 P(B) = 0.2041606337213415 P(O) = 0.2966569268290958 P(AB)= 0.09530792055695961 yes | ?- sampler(10,Samples). Samples = [a,o,o,a,b,a,ab,a,o,b] ?
(2) 隠れマルコフモデル
近年音声認識などの応用分野で広く用いられている隠れマルコフモデル(以下
HMM)は確率的オートマタであり、各状態において与えられた分布に従って記号
を出力する。外部からは出力記号のみ観測可能で、どのような状態遷移が行な
われたかは分からない。図 --6 に状態集合 {s0, s1}、
出力記号集合 {a, b} のHMMを掲げる。
長さ 5 の文を出力する HMMをPRISMプログラムで記述してみよう。時刻 Time, 状態 State において 記号 Obs を出力するという確率事 象 output( Time, State, Obs) と時刻 Time に 状態 State1 から状態 State2 に遷移するという確率事象 trans( Time, State1, State2) に対し、確率分布(各事象間 の統計的独立性を含む)を満たすように節のボディに bsw(i,n,x) を配する。統計モデル部分と制御宣言部分は 図 --7 のように記 述できる。ターゲットアトムは hmm(・) である。 hmm/1 は hmm3 を呼び出すが hmm/3 は再帰を引き起こす。このように PRISMは再帰呼び出しの記述を許すので我々は理論的には無限個の確率変数(確 率事象)を扱うことができる。
図 --6 の確率パラメータでサンプリングした500個の教師データを hmm.dat に与え学習させた。その結果表 --3 でグループID i のパラメータについて Conv の値に収束した(全ゴールの連言の対 数尤度の差が 10^-6 以下になったとき収束したと考える)。 Smp は サンプリングのときの確率パラメータである。
HMMには3つの基本問題とその解決アルゴリズムが存在する [1,11,18]。与えられた観測系列の確率を多 項式時間で計算する前向き確率計算アルゴリズム、観測系列が与えられたとき 状態の最尤パスを求める Viterbi アルゴリズム、HMMの訓練アルゴリズムであ る Baum-Welch アルゴリズムである。以下では前向き確率計算アルゴリズムを 実現するPRISMのユーティリティプログラムを記述してみよう。
○ 前向き確率計算アルゴリズム
観測列 O1O2…OT が与えられたとき、P(O1 O2…OT) を 計算する。まず
を用意し以下の計算を行なう。
ただし S は状態集合で S={s1,...,sN} である。また P(O|s) は状態 s のとき O を出力する確率(出力確率パラメータ)、P(s'→s) は状態 s から s'に遷移する確率(状態遷移確率パラメータ)を表す。 帰納部の計算の様子を 図 --8 に示す。
PRISM実行系の組み込み述語 prob/2 を使って
と確率は求められるが prob/2 は全解探索、すなわち全状態パスを求め 各パスの確率の和によって答を求めている。一方 図 --8 に明 らかなように前向き確率の計算アルゴリズムは αt (s) という変数に 時刻 t までの計算結果を格納しそれ以前の再計算は行なわず、効率的なア ルゴリズムになっている。PRISMに組み込みの prob/2 確率計算ユーティ リティは汎用性はあるが、その計算量/計算時間においてHMMのようなモデル とそれに特化(specialized)されたアルゴリズムの方が計算の効率は優れてい る。| ?- prob(hmm([a,a,b,a,b]),P). P = 0.054776591623507186 ?
このような計算効率に関する問題は今後の課題であるが、現段階においても PRISMの記述力により特化されたアルゴリズムをユーティリティプログラムと して記述することは可能である。次にその例を示す。
○ ユーティリティプログラム
図 --6 のHMMに対する前向き確率の計算アルゴリズムを図 --9 にユーティリティプログラムとして記述しよう。すると 次のように前向き確率が計算される。
図 --9 中には実行系の組み込み述語 prob/2 が使われ ている。ただしこのプログラム中では output/3, trans/3 のよ うに bsw/3 を直接呼び出す述語の確率を求めるために使われているので、 探索の手間はとらない。 sum( List, Sum) はリスト List 中の全要素の和が Sum となることを表し、 member( Element, List) は Element が List にあることを示す。 findall( Templete, Goal, List) は List が Goal を満たすような Templete の代入例からなるリストである ことを示す| ?- forward([a,a,b,a,b],P). P = 0.054776591623507186 ?
ここでは省略するが Viterbi アルゴリズム、Baum-Welch アルゴリズムのユー
ティリティプログラムも同様にプログラミングできる。組み込みのEM学習ルー
チンの代わりに Baum-Welch アルゴリズムを使うこともできる。
(3) ベイジアンネットワーク
ベイジアンネットワーク(Bayesian network)は命題間の確率的因果関係を有向
非循環グラフ(DAG)で表現したものである。他に Belief network, Causal
network, Influence network などとも呼ばれている。ネットワークの形状か
ら各命題間の統計的独立性が規定され、独立性から全命題の同時分布が条件つ
き確率の積に還元される。例えば 図 --10 のベイジアンネットワー
クに関して次が成り立つ。
図 --10 のベイジアンネットワークを Poole [17] の 例に習いPRISMプログラムで記述する。ここでは観測者は報告( Report) と煙( Smoke)しか観測できないとし、ターゲット述語は world/2 に指定している。PRISMプログラムは
c_smoke(yes,Fi):- bsw(sm(Fi),none,1).のように bsw/3 の第1引数に変数 Fi を用いることを許している
% Control declarations: target(world,2). data('world.dat'). % Only Smoke and Report are observable. world(Sm,Re) :- world(_,_,_,Sm,_,Re). % Joint distribution: world(Ta,Fi,Al,Sm,Le,Re) :- fire(Fi),tampering(Ta),c_smoke(Sm,Fi), c_alarm(Al,Fi,Ta),c_leaving(Le,Al),c_report(Re,Le). % Marginal probabilities: smoke(Sm) :- fire(Fi),c_smoke(Sm,Fi). alarm(Al) :- fire(Fi),tampering(Ta),c_alarm(Al,Fi,Ta). leaving(Le) :- alarm(Al),c_leaving(Le,Al). report(Re) :- leaving(Le),c_report(Re,Le). % Conditional probabilities: tampering(yes):- bsw(ta,none,1). tampering(no) :- bsw(ta,none,0). fire(yes):- bsw(fi,none,1). fire(no) :- bsw(fi,none,0). c_smoke(yes,Fi):- bsw(sm(Fi),none,1). c_smoke(no,Fi) :- bsw(sm(Fi),none,0). c_alarm(yes,Fi,Ta):- bsw(al(Fi,Ta),none,1). c_alarm(no,Fi,Ta) :- bsw(al(Fi,Ta),none,0). c_leaving(yes,Al):- bsw(le(Al),none,1). c_leaving(no,Al) :- bsw(le(Al),none,0). c_report(yes,Le):- bsw(re(Le),none,1). c_report(no,Le) :- bsw(re(Le),none,0).
そしてサンプル値 Smp を使って1000個サンプリングし、これを教師デー タとして与え、学習させた。その結果452回の繰り返しで収束し、表 --4 の結果を得た。表 --4 を見ると個々の収束 値 Conv はサンプル値 Smp と異なっているが、観測可能なレベ ルでの確率分布 PDB(world(・,・)=1) は十分近付 いている。不十分な情報を補い推定するEMアルゴリズムの特質がよく現れてい るといえる。
(4) 回路故障モデル
Poole の確率的 Horn アブダクション言語 [17] の例として引用さ
れていた de Kleer の3つのインバータの直列回路 [3]の故障モ
デルの例を掲げる(図 --11)。各インバータは良好( ok)、
ショート( shorted)、断線( blown)の3つの内部状態をもち、これ
らの内部状態はある確率分布に従っているとする。我々は再びPoole に習い、
この回路の故障モデルをPRISMプログラムで記述する。ターゲット述語は
circuit/2 で、 circuit( In, Out) とすると、外部からは
回路の入力 In と出力 Out のみが観測可能である。
インバータの内部状態を良好、ショート、断線がそれぞれ 0.8, 0.05, 0.15 という確率分布、かつ入力 on/ off は一様分布に従って(ただし 内部状態と入力は独立事象とする)サンプリングしたものを教師データとして three_inv.dat に格納する。学習を行なうと1005回の繰り返しで表 --5 の収束値 Conv を得た。また表 --5 より ok(g), shorted(g), blown(g) が真になる確率は十分近付いたと考えられる。EMアルゴリズムに より外から観測できない g1-- g2 間、 g2-- g3 間の 値を補いながら推定が行なわれている。
circuit(In,Out) :- val(in(g1),In),val(out(g3),Out). val(out(G),on) :- ok(G),val(in(G),off). val(out(G),off) :- ok(G),val(in(G),on). val(out(G),V) :- shorted(G),val(in(G),V). val(out(G),off) :- blown(G). val(in(G),V) :- conn(G1,G),val(out(G1),V). conn(g1,g2). conn(g2,g3). ok(G) :- bsw(s1,G,1). shorted(G) :- bsw(s1,G,0),bsw(s2,G,1). blown(G) :- bsw(s1,G,0),bsw(s2,G,0). target(circuit,3). data('three_inv.dat').
本研究では以上に述べたようなPRISM処理系の試作版を SICStus Prolog 2.1 (Cルーチン呼び出し使用)で実装した。ソースプログラムは2000行程度である。 現在ユーザインタフェースは Prolog に委ねている。また使用計算機/使用OS はそれぞれ Sun SparcStation 10/SunOS 4.1.4 を想定している。ドキュメン トとして本報告書で述べたようなPRISMに関する概論、組み込み述語の使用法 およびプログラム例を含むプログラミングガイドからなるユーザ用マニュアル を添付する。
Prolog-like の形で統計モデルの記述とパラメータ学習が可能になったので、 研究としての成果として述べたことはそのままソフトウェアとしての成 果となるだろう。加えてユーザはユーティリティプログラムとして統計モデル をメタ的に扱うことができるので、前節の例のようにモデル中に付与されたパ ラメータの学習を行なうだけでなく、例えばモデル構造を決定するようなプロ グラムを書き、決定を司るアトム(もしくは節)に bsw/3 を配すること で構造を学習するプログラムも書けるだろう。
前述のようにユーザはPRISMプログラムを Prolog に幾つかの組み込み述語が
加わった程度の意識で書けることが見込まれるが、加えてPRISM処理系はイン
タフェースをユーザにとって使いやすいものになっている。PRISM処理系がユー
ザに有益な機能を分かりやすい形で提供することを以下のPRISM処理系の動作
の様子をもって示そう。
(1) 実行系
実行系の動作例として先に挙げた血液型遺伝子モデルをPRISM処理系に読み込
み、確率計算のツールとして活用してみよう。各 bsw(i,n,x) パ
ラメータの値は既に学習(もしくは人工的に設定)されているものとする。
という質問するとターゲットアトム bloodtype(・) の確率分布に 従って| ?- sample(bloodtype(X)).
X = a?
, X = b?
, X = o?
,
X = ab?
のいずれかの答を返す。また
という形の質問であれば実行系はターゲットアトムの確率分布に従って yes または no を返す。| ?- sample(bloodtype(ab)).
また gene(dad,a) が真の時に bloodtype(a) が真である条件つ き確率は| ?- prob(bloodtype(a),P). P = 0.38971075015247592 ?
より各々求められる。また prob/1 は prob( Formula) と して与えられたとき答えの確率値を変数に束縛する代りに Formula を 清書して出力する。例えば| ?- cprob(bloodtype(a),gene(dad,a),CP). CP = 0.83697701434085235 ?
また prob( Formula) において式 Formula を真にする(説 明する)ような bsw(・,・,・) の連言が見つからない ときは次を出力する。| ?- prob((bloodtype(a);bloodtype(b))). The probability of bloodtype(a) v bloodtype(b) is 0.6080351526139445. yes
| ?- prob((bloodtype(a),bloodtype(b))). bloodtype(a) & bloodtype(b) is false with probability 1. yes
各々の bsw(i,n,x) は独立であるため確率を計算するには論理積 で確率積、論理和で確率和を計算すればよい。
例えば bloodtype(a) に対する確率計算式は次のようになる。 v は論理和 ∨、 & は論理積 ∧ を表す。
probf/2 では probf( Formula, DNF) としたとき DNF に bsw(i,n,x) の選言標準形を表すリストを 返す。| ?- probf(bloodtype(a)). bloodtype(a) is explained by bsw(1,dad,1) & bsw(1,mum,1) v bsw(1,dad,1) & bsw(1,mum,0) & bsw(2,mum,0) v bsw(1,dad,0) & bsw(1,mum,1) & bsw(2,dad,0)
| ?- probf(bloodtype(a),Formula). Formula = [[bsw(1,dad,1),bsw(1,mum,1)], [bsw(1,dad,1),bsw(1,mum,0),bsw(2,mum,0)], [bsw(1,dad,0),bsw(1,mum,1),bsw(2,dad,0)]] ?
(2) 学習系
学習系の述語は learn/0 または learn/1 である。節
(1) 血液型遺伝子モデル の血液型遺伝子プログラムの学習の様子を以下
に示す。学習系の処理の様子は逐一ユーザに知らされる。
| ?- learn. {Extract goals from blood.dat} ← blood.dat から教師データを取り込む {building explanations...} ← 対応表を構築 .......................................................................... .......................... {built explanations for all the goals.} BSW 1: 0.752654 BSW 2: 0.449459 ← 各 bsw パラメータをランダムに設定 [0] Log_like= -210.984965 ← 全ゴールの連言の対数尤度 BSW 1: 0.396859 BSW 2: 0.357359 [1] Log_like= -136.017496 ( パラメータ値の更新を繰り返す ) : BSW 1: 0.292315 BSW 2: 0.230360 ← bsw 1 のパラメータ値 0.292315 に収束 [11] Log_like= -128.004797 bsw 2 のパラメータ値 0.230360 に収束 -- converged! ← 終了 yes
(3) その他
PRISM処理系はプログラミングを支援する目的で実行/学習系の組み込み述語 の他にもパラメータ操作、サンプリングユーティリティなどに関して組み込み 述語を用意している。
一方 set_bsw/1, set_bsw/2 は bsw(i,・,・) のパラメータ θi の値を設定するもので ある。学習系では教師データからパラメータ値が決定するが、 set_bsw/1, set_bsw/2 を用いれば人手で設定できる。 set_bsw/2 では個々のパラメータ θi を設定し、 set_bsw/1 では記号 : を使って複数のパラメータを同時に設定する。すなわち次 の2つは同様の効果をもつ。| ?- show_bsw.
save_bsw/1 は save_bsw( Filename) の形でファイル名 Filename のファイルに bsw(i,・,・) のパラメー タ値 θi をセーブする。逆に restore_bsw/1 は restore_bsw( Filename) の形でファイル Filename からパラメー タ値を復帰させるものである。| ?- set_bsw(1,0.6),set_bsw(2,0.3). | ?- set_bsw([1:0.6,2:0.3]).
とすれば X = head ?, X = tail ? のいずれかの答が等確率 (0.5)で得られる。| ?- dice_uniform([head,tail],X).
dice_multi/3 は dice_uniform/2 と似ており、 dice_multi( Atoms, Probs, Result) の形で定数項のリス ト Atoms とそれが出現する確率のリスト Probs (その要素の和 は 1 でなければならない)が与えられたとき、その確率分布に従って Atoms の中から1つ選び、 Result に束縛する。例えば
とすれば X = a ?, X = b ?, X = o ?, X = ab ? という答がそれぞれ確率 0.4, 0.2, 0.3, 0.1 で返ってくる。 dice_uniform/2, dice_multi/3 を用いたプログラミング手法 を次に示す。| ?- dice_multi([a,b,o,ab],[0.4,0.2,0.3,0.1],X).
まず制御宣言で次のように教師データファイルを特別なファイル名
user にする。すると翻訳系により learn/1 が生成される。
learn( Goals) は Goals に教師データのリスト(たとえば
[bloodtype(a),bloodtype(b),...]
)を受けとり学習を行なうコマンド
である。
次に以下をユーティリティプログラムとして加える。プログラム中の japanese_sampler/2 は japanese_sampler( N, Goals) の形で N が与えられたときに N 個の bloodtype(・) のサンプリングデータを取り出すものである。data(user).
my_learn(N) :- japanese_sampler(N,Goals), learn(Goals). japanese_sampler(0,[]). japanese_sampler(N,[bloodtype(X)|Goals]) :- N > 0, dice_multi([a:0.4,b:0.2,o:0.3,ab:0.1],X), N1 is N-1, japanese_sampler(N1,Goals).
書き換えた血液型遺伝子プログラムをPRISM翻訳系に読み込ませた後、次のよ うに命令すると、100個のサンプリングデータを教師データとしてパラメータ 学習を開始する。
| ?- my_learn(100).
これまで見てきたようにPRISMプログラミング処理系は全く新しいタイプの論 理プログラミング処理系である。従って課題は山積みになっているといってよいだろ う。例えば以下のようなものが考えられるが、理論的整備と実装を含め全てを 実現するにはかなりの時間が必要だろう。2年目の研究ではこの内の幾つかを 選んで実現し、ソフトウェアとしてのPRISMに反映させたい。
まずPRISMプログラムのクラスの拡張があげられる。現在BSプログラムでは i≠j または n≠m のとき bsw(i,n,・) と bsw(j,m,・) は互いに独立な確率的スイッチだった。それに対しこ の2つのスイッチが相関をもつような確率分布、例えば Boltzmann 分布をもつ ようなBSプログラムが考えられる。また x に 0, 1 の2値しか 許さない bsw(i,n,x)に対し多値( multiple values)を許す 確率的スイッチ msw(i,n,x) をもつMSプログラムが考えられる。 それぞれ学習アルゴリズムなどの理論的見直しを経てPRISM上への実装を目指 したい。
現在のPRISM言語の設計では1つのターゲット、1つの教師データに対して大局 的(global)学習を行なっている。しかし局所的(local)な学習を行なってそれ を組織的に束ねることができればまた違った応用の見通しが出てくるだろう。 このような局所的な学習を考慮するようにPRISM言語を設計し直す必要がある。
また現在のPRISMの学習方式は事前にある程度の大きさの教師データを用意し、 そのデータに振舞いを合わせるというバッチ学習を採用していた。それに対し オンライン学習と組み合わせるというアプローチが考えられる。
実用的な計算時間を考えた場合、学習系における全解探索の計算時間への負荷 が問題になる。サンプリング探索、並列探索などを用いた全解探索に代わる効 率的な探索手法を模索する必要がある。
最後にユーザとの関わりで見た場合以下のようなプログラミング環境の整備が 必要になる。これらは2年目の研究で行なう予定である。