組み合わせ数をグラフパスとして数え上げる

[statgenetJimu]

ZDDではパスの数え上げができますが、それを使って「正確確率検定」をZDD的に考えてみることにします。たとえば n x m の長方形格子の左下隅から右上隅へのパスの数は、右に行くか上に行くかのみを許す(directed graphとみなす)ならば、『全部でn + m歩あるくパスであって、そのうち、右にn歩、上にm歩であるパス』のことなので、結局、(n+m) C n となることになります。 2 x 2 分割表でそうサンプル数Nの周辺度数が　n+m=N, a+b=Nとなっているときには、常に４つの移動が許されたグラフでのパスを想定すればよいです。４つの移動をx,y,z,wとすると、総歩数 N = x+y+z+wであり、周辺度数制約からx+y=n,z+w=m,x+z=a,y+w=bという制約のあるパスを考えればよいです。 それは、x+y=n,z+w=m,x+z=a,y+w=bを満足する複数のゴール地点へのパスの数を列挙することになります。

もちろん、こんな面倒なことをしなくても組み合わせ数計算でできるわけですが、「分割表に様々な条件が入ったりしたときにもグラフとその上のパス」とに話が落とせれば、正確確率が計算できることになりそうです。

さて。Graphillionでは無向グラフを対象にしているようですが、どうしたら、n x m の長方形エリアの左下隅から右上隅へのパス数をZDD的に扱えるでしょうか。

[skoyama427]

p03.py

n*mの格子を左下から右上までn+m歩で歩くということは、この間を最短距離で歩くということに等しいと思います。これはすべてのグラフ集合のうちエッジの数がn+mのものを抜き出すということと同じです。graphsetオブジェクトのクラス関数にエッジ数で指定したグラフを取り出すgs.graph_size()というものがあるのでこれを使いました。

from graphillion import GraphSet import graphillion.tutorial as tl n=3 m=2 nm = n+m goal = (n+1)*(m+1)

universe=tl.grid(n,m) GraphSet.set_universe(universe)

Gt=GraphSet.paths(1,goal)

Gs=Gt.graph_size(nm) example = Gs.choice() example = (Gs.including(5)).choice()

print len(Gs) print example tl.draw(example)

[statgenetJimu]

ああ、そうですね。 p03.pyは/python/src/に上げておきました。 じゃあ、n x m格子グラフがあったときに、それぞれのエッジに値を持たせるとしましょう。 その上で、列挙されたパスについて、パスに含まれるエッジの値の和をとることはできるでしょうか。 もしそれができるとすると、各パスに「そのパスを取る確率の対数」を与えることで、それぞれのパスが発生する確率の対数が計算できそうです。

その上で、次のような例を考えます。 今、治療法Aと治療法Bとがあり、それぞれ成功率がPa,Pbとします。N人の患者をランダムにA、Bに割りつけながらN人目まで到達するとき、「１人目はAで成功、２人目はAで失敗、３人目はBで失敗・・・」というようなパスができますが、その生起確率が計算できるでしょう。 ランダムアサインメントだとつまらないのですが、まずはそれをやってみるとして、次のように考えます。k人目の人にA、Bのどちらを採用するかを、k-1人目までの成否の人数情報（A,Bの２治療法の２つの帰結があるので４つの整数）で決めることにします。 簡単には、As、Af、Bs、Bfの４整数で、As/(As+Af) > Bs/(Bs+Bf) ならAをとり、その逆ならBを取り、等しければ0.5の確率でAを取る、というようなルールも入れられそうです。このA,B選択ルールに工夫を加えると面白いのですが。。。 そんなことができるでしょうか。 実際、これは、希少疾患等での決断理論用として、院生のWangさんがこの確率の性状についての研究をしているのですが、現在は、ZDDを使わずに計算しています。

[statgenetJimu]

start -> goal の部分グラフを列挙してから、その長さで選択をかけるより、指定の長さの部分グラフを列挙してから start->goalの条件で選択するほうが速い（最初に作る部分グラフ集合の要素数が少ない？もしくはZDD化した構造が簡単だから？）ようなのでp03RY.pyとしてアップしておきました。

[skoyama427]

p04.py Eki.py内のsumWeightメソッドなどをみると、パスの重みはそのまま扱えるような形では組み込まれていないようです。簡単ですがp04.pyに重み付きの格子を作るメソッドを定義してみました。

p03RY_2.py 最初にできる部分グラフ集合はstart -> goalで選択したもののほうがsize=mnよりも小さいようです。にもかかわらず実行時間が大きいのはZDDの複雑性が違うのでしょうか。ノードIDなどで複雑さが想定できそうですのでまた調べてみます。

[statgenetJimu]

重み付き計算については、もう少し考えてみたいと思います。 GraphSet.pathsとGraphSet.graph_sizeとの速さの違いは

pathsの方が return GraphSet.graphs(vertex_groups=[[terminal1, terminal2]], degree_constraints=dc, no_loop=True, graphset=graphset) と graphsという関数を使っているのに対し、 graph_sizeの方が return GraphSet(self._ss.set_size(size)) とpythonに作り込まれている集合クラスを使っているのが大きいのかもしれません。

---

単なるnCmとかの数え上げだと、重いばかりで何の嬉しいことも思いつかないので、なにか複雑な条件付き数え上げ→正確確率→検定、というシチュエーションがないか、検討したいと思います。