index.html

Combinatorica でネットワーク解析を学ぼう

前口上

この文書は吉野がネットワーク解析について学ぶために作成した私的なメモである．従って，このメモの誤りによって生じたいかなる不都合もその責任はメモに頼った人が責任を負うべきものであり，吉野は一切の責任を持たない．しかし，このメモを更に発展させてほしいという要望は作者のメモ作成の駆動力となるので，誤りの報告・追加の要望・暖かい激励の言葉は takashiあっとまーくrandom-walk.org に送信してみる価値はあるかもしれない（「あっとまーく」は半角の＠で）．

このメモを書くきっかけになったのは Watts and Strongatz と Albert and Barabasi である（参考文献参照）．これらの論文の内容を自分なりに確認することと Combinatorica について学ぶことがこのメモの主題である．そのため，一般のグラフ理論の本よりもかなり偏った内容になっていると思われる．そのあたりは個人的な趣味の話なので勘弁してほしい（汗）それにずっと作りかけのままだし．次の目標はランダムグラフの相転移あたりかもしれない．

この文書の更新履歴は以下のとおり．
・2003/12/10 Ver. 1.0
・2004/02/20 Ver. 1.1 (パーコレーションを追記）
・2004/08/21 Ver. 1.2 (ゲーム理論の部分を追記，Small world 部分を変更）
・2004/10/06 Ver. 1.3 (複数の場所でミスを発見，small world 部分に正方格子の場合の話題を追加)
・2005/01/10 Ver. 1.4 (リクエストにより Scale Free Network のプログラムを追加，説明文のミスを修正)
・2005/01/24 Ver. 1.5 (パーコレーション部分を追加・修正）
・2006/10/19 Ver. 1.6（リクエストにより Scale Free Network のグラフィックスに正多角形バージョンを追加．説明文のミスを修正．スタイルを修正）
久しぶりに眺めてみました．今となっては恥ずかしい幾つかのミスや表記があります．修正できるところは修正しましたが，今の自分ならこうは書かないなという部分も多いです．コメントをいただけると修正できてうれしいです．単純なことでもいいので何かありましたらご指摘願います．（2006/10/19 追記）

テクニック

基本的なこと

パッケ-ジのロード

はじめに Combinatorica パッケージをロードしなくてはならない.

とする．これは,毎回の作業の始まりに一度だけ行えばよい手続きである．頻繁に使うようになったら， Mathematica を起動したときに自動的に読み込むようにしておくとよいかもしれない．

Graph オブジェクト

Combinatorica でのグラフの内部表現は変えない方が，用意されているさまざまな関数を使いまわすことができるので好都合である．それがどのような要素で構成されているのかを確認するために，まず， 5つの点からなる完全グラフを作ることにする．完全グラフというのが何かは後で考えることにしよう．

このように,残念ながら作られたグラフがどのようなものなのかは,その表示を見るだけでは本当の概要（点と辺の数そして有向・無向）しかわからない. FullForm を用いてすべて表示すると，このオブジェクトはList の組み合わせからなることがわかる．

グラフを表現するために Combinatorica が保持しているデータは連結に関するデ-タと位置に関するデータであることがわかる．ヘルプを参照すると，

と，辺と点座標のリスト＋オプションであることがわかる．

グラフの表示

上記のように，グラフオブジェクトを作成するだけでは，その構成を確認するのは難しい．グラフを表示するには，

[Graphics:HTMLFiles/index_11.gif]

とする． ShowGraph には多くのオプションが用意されている．これを確認する．

例えば点の色を変えたい場合には,

[Graphics:HTMLFiles/index_15.gif]

さらに個々の点や線のスタイルのみを変えることもできるのだが,これは本題から外れるので Combinatorica のヘルプや関連サイトで情報収集をしてほしい.

主な解析ツ-ル

点(vertex)の数を数えるには V を用いる．たとえばの頂点の数は,

で求められる.同様にして，辺（edge）の個数を求めるには， M を用いる．

辺のリストは Edges で得ることができる．

隣接の関係（点と点が繋がっているのかいないのか）は，隣接行列を用いて表現される．

( {{0, 1, 1, 1, 1}, {1, 0, 1, 1, 1}, {1, 1, 0, 1, 1}, {1, 1, 1, 0, 1}, {1, 1, 1, 1, 0}} )

数字は繋がっている辺の本数を表している．このグラフは向きを持っていないグラフ（無向グラフ）なので，転置行列も同じ形をしている．対角成分がすべてゼロであることは自分自身へのループが存在してないことを意味している.

頂点の個数は次数と呼ばれる．次数の平均を位数と呼ぶことがある．これは辺の数の２倍を点の数で割ればいいので，

と求められる．個々の点についての位数を求める場合には，隣接行列を用いる．隣接行列は List の List であることを確かめておく．

これを用いると，それぞれのリストは各点から出ている辺の個数に等しいので，各リストに Plus を適用することによって，各点の位数が得られる．

2点間の経路は ShortestPath 関数で求められる．例えば k5 において，点２と点５間の経路を求めたい場合には，

その距離は経路のリストの長さから1を引けばいいので,

で求められる.

すべてのペアについて距離を求めたい場合には， AllPairsShortestPath を用いる（完全グラフの場合には，２点間の距離はすべて 1 である）．

( {{0, 1, 1, 1, 1}, {1, 0, 1, 1, 1}, {1, 1, 0, 1, 1}, {1, 1, 1, 0, 1}, {1, 1, 1, 1, 0}} )

グラフの種類

完全グラフ

[Graphics:HTMLFiles/index_39.gif]

である.

木と星

木（Tree）とは閉じたループも孤立した点もひとつも持たないグラフをいう.そのためには (n-1) 本の線で構成されていなければならない．ただし，その方法は一意的でないために， Combinatorica は Tree という関数を持たない．線の配置をランダムに指定しても良い場合には， RandomTree を用いる．たとえば，

[Graphics:HTMLFiles/index_41.gif]

という要領である.星(Star)とは木の中で，ひとつの点が (n-1) の線を持ち，すなわち位数が(n-1)であり，他の点がひとつしか線を持たない，すなわち位数が 1 であるグラフをいう．これはグラフが一意的に決定されるので， Combinatorica で関数が用意されている．

[Graphics:HTMLFiles/index_43.gif]

平均距離はほとんどが距離２にあるので，

となる．当然のことながら頂点の数が大きくなるほど２に近づいていき，

[Graphics:HTMLFiles/index_47.gif]

規則的な木を作りたい場合には CompleteKaryTree を使うとよい．

$ShowGraphArray[Table[CompleteKaryTree[(2^i - 1)/1, 2], {i, 2, 5}]] ;$

[Graphics:HTMLFiles/index_50.gif]

$ShowGraphArray[Table[CompleteKaryTree[(4^i - 1)/3, 4], {i, 2, 4}]] ;$

[Graphics:HTMLFiles/index_52.gif]

$ShowGraphArray[Table[CompleteKaryTree[(8^i - 1)/7, 8], {i, 2, 4}]] ;$

[Graphics:HTMLFiles/index_54.gif]

$Table[(Plus @@ (AllPairsShortestPath[CompleteKaryTree[2^i - 1, 2]]//Flatten) )/(i * (i - 1)), {i, 2, 8}]//N$

サイクル

サイクルとはひとつの閉じたループのみで構成されるグラフである.

[Graphics:HTMLFiles/index_58.gif]

Wheel

星とサイクルをあわせて車輪（Wheel）が作られる．

[Graphics:HTMLFiles/index_60.gif]

Empty Graph

Empty graph とはひとつも線を持たないグラフを言う.

[Graphics:HTMLFiles/index_62.gif]

グリッドグラフ

視覚的に言えば格子状に配置したグラフである.でも,この言い方はトポロジカルな関係だけを問題にするグラフ理論では適切な言い方ではないかもしれない．でも，直感的には一番納得しやすいのではないかと思う．２次元や３次元の正方格子に相当するグリッドグラフはあるみたいだが，多次元の正方格子に相当するグリッドグラフはないのだろうか．作ろうと思えば作れるのだが…あるとパーコレーション理論との関係を検証できて便利なんだけど…

[Graphics:HTMLFiles/index_64.gif]

ランダムグラフ

ランダムグラフは,このノ-トを書く動機になったグラフである.これはn個の点からなり,そこに含まれる任意の2点をつなぐ線の有無を確率 p で決定することによって作られる.このグラフの持つ性質については次の章でまとめるつもりである．Combinatorica では RandomGraph とExactRandomGraph というふたつの関数が準備されている．前者は確率を後者は陵の本数を与えてグラフを作成する．すなわち，RandomGraph[n, p] は n 個の節点によって構成される完全グラフをのすべての陵を確率 p で生き残らせるか消去するかを決定することによって与える．それに対して ExactRandomGraph[n,e] はe本の生き残る辺をランダムに選択する．RandomGraph にはさらに向きの有無なども指定できる．

[Graphics:HTMLFiles/index_67.gif]

[Graphics:HTMLFiles/index_69.gif]

RandomGraph は与えられた確率に従ってグラフを作るだけなので，作られた稜の数が一定ではない．

Circulant グラフ

Circulant グラフは周期的に連結するグラフを作る．Circulant[n,1]は点の数が n で隣と連結するグラフ，すなわち Cycle グラフである．

[Graphics:HTMLFiles/index_73.gif]

Circulant[n,k]は k だけ飛んだ先にある点と周期的に連結するグラフを作る．

[Graphics:HTMLFiles/index_75.gif]

正則グラフ

すべての節点の位数が等しいグラフを正則グラフという．位数が k である正則グラフをk次の正則グラフという． Combinatorica には対称性の少ない連結，すなわち連結の仕方がランダム（ヘルプの言葉を借りれば「準ランダム」）な正則グラフを作る関数 RegularGraph が用意されている．

[Graphics:HTMLFiles/index_78.gif]

これは呼び出されるたびに異なる正則グラフを作る．最も対称性のよい正則グラフを作るには，CirculantGraphを重ね合わせればよいので，以下の関数を使う．

[Graphics:HTMLFiles/index_82.gif]

自分でグラフを作るには

既成のグラフだけで足りることばかりではないと思われる．そこで，点や線を付けたり消したりする方法を知っておくと便利である.

グラフ g に新しい辺を加えるには AddEdges[g, edgeList] を用いる．ここで edgeList は連結するふたつの点のリストである．例えば，前述のランダムツリーに線をひとつ加えてみよう．初めに，その構成を確認しておこう．隣接行列を使って隣接の様子を再確認する．

この行列の要素のうち値がゼロであるものが隣接していないことを意味している．例えば節点 4と接点 8，接点3と接点20がそれである．

[Graphics:HTMLFiles/index_89.gif]

が隣接させる前のグラフである．AddEdge を用いて，２点を隣接させてみる．

nrt20 = AddEdges[rt20, {{4, 8}, {3, 20}}] ; ShowGraph[SetGraphOptions[nrt20, {{4, 8 , 3, 20, VertexColor→Green}, {{4, 8}, {3, 20}, EdgeColor→Red}}], VertexNumber→On] ;

[Graphics:HTMLFiles/index_91.gif]

逆に削除することもできる．

ShowGraph[SetGraphOptions[nnrt20, {{4, 8 , 3, 20, VertexColor→Green}, {{4, 8}, {3, 20}, EdgeColor→Red}}], VertexNumber→On] ;

[Graphics:HTMLFiles/index_94.gif]

同様にして点の追加も行える.

AddVertex[nnrt20]//ShowGraph[SetGraphOptions[#, {{4, 8, 3, 20 , VertexColor→Green}, {{4, 8}, {3, 20}, EdgeColor→Red}}], VertexNumber→On] & ;

[Graphics:HTMLFiles/index_96.gif]

点を削除することも可能である.

DeleteVertex[nnrt20, 20]//ShowGraph[SetGraphOptions[#, {{4, 8 , 3, 20, VertexColor→Green}, {{4, 8}, {3, 20}, EdgeColor→Red}}], VertexNumber→On] & ;

[Graphics:HTMLFiles/index_98.gif]

大まかな形を作っておいて,そこから副グラフを作ることで新たなグラフを作ることもできる.

ふたつ以上ののグラフの和集合を作る．

[Graphics:HTMLFiles/index_100.gif]

あるグラフから一部（副グラフだ）を取り除く．

[Graphics:HTMLFiles/index_102.gif]

副グラフを作る．副グラフを作ることは生き残りを選ぶことといえるかもしれない．もとのグラフと生き残る点のリストを与えて，その他の点とその点に繋がっていた辺を除くことができる．

[Graphics:HTMLFiles/index_105.gif]

[Graphics:HTMLFiles/index_107.gif]

解析の方法

グラフ（またはネットワーク）を特徴付ける量には上に述べたような点や辺の数や位数などの他にもさまざまな量がある．その中でもネットワークの構造解析を行うために使われることが多い，直径（diameter)や平均最短経路，そしてクラスター係数（clustering coefficient）の Mathematica による計算方法について説明する．

直径

以下の文は Diameter が実装されていたことをを知らないときにかいたもの．

直径は既存の辺を通って２点間を連結する最小値（最短経路）の最大値である．最短経路は， AllPairsShortestPathを用いて計算されるので，その最大値を求めればよい．
ネットワークの構造解析という観点から言うと，直径が小さいことは遠くに行くまでにかかる手数が少なくて済むことを示しているので，小さいほどそのネットワークが密であると言うことができる．

[Graphics:HTMLFiles/index_109.gif]

（Max は Flatten 属性を持っているから Flatten は省けるか）以下のような関数を作っておけば便利かもしれない（これも Flatten が省けるな）．

ランダムグラフのリストを作って，直径の平均値の確率依存性を調べてみる．

[Graphics:HTMLFiles/index_125.gif]

以下はDiameter を知った後のコメント．
おそらくDiameter の実装は上の実装とほとんど変わらないのではないかと思われる．大きなネットワークだと時間がかかりすぎるのが難点だ．

平均最短経路

直径と似たような性質を議論できる量として，平均最短経路がある．これは２点間の最短経路の平均である．すなわち，AllPairsShortestPath の要素を２点を選ぶ組み合わせ数の２倍で割った値である．

従って以下の関数で表現できる．

これも，ランダムグラフのリストを作って，平均値の確率依存性を調べてみる．

[Graphics:HTMLFiles/index_132.gif]

クラスター係数

クラスター係数（clustering coefficient）は周囲の点との関係の深さを表す指標である．位数 k の点は連結しているすべての点が互いに連結しているとき， k(k-1)/2個の辺で連結されることになる．この k 個の点のうち，実際に連結している辺の数の k(k-1)/2 に対する辺をクラスター係数と呼ぶ．この値を隣接行列から求める方法を考える．

ひとつの行について１がある列をピックアップしてみる．

この集合から作られるペアは,

なので，これらが連結している（すなわち要素が１）なのか，連結していない（すなわち要素が０）なのかを確かめればよい．これは各要素を抽出するだけの話なので，

これらの和を k(k-1)/2 で割った値が求めるクラスター係数である．

ここまでの話を整理して，すべての点について一度にクラスター係数を求められるような関数を作成してみることにする．

まず，各点における位数を求めるには，前述のように，隣接行列の各行について和をとればよいので，

で良い．これらを k(k-1)/2 倍するには，

でよい．

次に各点についてクラスター係数の分子を求める方法を考える．上に述べた方法を順次あてはめていくには，

counts = Apply[Plus, Map[Extract[adRg2004, #] &, Map[KSubsets[#, 2] &, Map[Flatten, Map[Position[#, 1] &, adRg2004], {1}], {1}], {1}], {1}]

このようにして作られたふたつのリストを互いに割ればいいので，

と求められる.その平均は,

以上の結果もひとつの関数にまとめられる．ただし位数がゼロの場合を考慮してないことに注意をする．

こいつも，ランダムグラフのリストを作って，平均値の確率依存性を調べてみる．

[Graphics:HTMLFiles/index_156.gif]

ランダムグラフ

ここではランダムグラフの性質について更に詳しく書く予定．続きを読みたいと思ったらリクエストメールを送ろう（笑）

[Graphics:HTMLFiles/index_159.gif]

ランダムグラフの集合を作って,その統計的な性質を見てみる.

辺の本数を数えて平均を取ってみる.

解析解は以下のとおり.

このグラフの集合を用いて幾つかの性質を確認してみよう.例えば，ひとつのグラフについてその直径を求めるためには，Path の中で最も大きい値を求めればよい．

Small World Network

最近流行りの Small World Network の性質を Mathematica を使って調べてみよう.以下ではグラフと同じ意味でネットワークという言葉を用いる．

Small world の構成方法

Small World と呼ばれるグラフには何種類かある．ここでは Watts と Strongatz の論文に書かれた方法とほぼ同じ方法で Small World Network を構築することを試みる．

はじめに規則的なネットワークを作成しよう．前述の正則グラフを用いる．

[Graphics:HTMLFiles/index_174.gif]

ひとつの辺を指定して，それを現在連結されてない二点を結ぶ辺に置き換える関数を作成する．

[Graphics:HTMLFiles/index_178.gif]

[Graphics:HTMLFiles/index_179.gif]

[Graphics:HTMLFiles/index_180.gif]

更にこれを全ての辺について適用するために Fold を用いる．

[Graphics:HTMLFiles/index_182.gif]

Fold する前に変更の有無を選択する関数を作る．

これを用いれば Small World の完成．

[Graphics:HTMLFiles/index_186.gif]

ひとつの関数にまとめておく.

[Graphics:HTMLFiles/index_190.gif]

[Graphics:HTMLFiles/index_191.gif]

[Graphics:HTMLFiles/index_192.gif]

[Graphics:HTMLFiles/index_193.gif]

[Graphics:HTMLFiles/index_194.gif]

p が大きくなるとともに半径（の平均値）が小さくなることがわかる．

以下， Small World の性質についても調べたいが，これも次の更新に譲ることにしよう．リクエスト待ってます（来ないってｗ）

グリッドグラフに small world を適用した場合

周期境界条件を課した二次元正方格子上での small world を作成する．

はじめにグリッドグラフを作る．

[Graphics:HTMLFiles/index_202.gif]

周期境界条件を用いたいので， n × n のグリッドグラフに対して， i n 番目と ( in -4) 番目，および 1 から n 番目と n(n-1) +1 から n(n-1) + n までをそれぞれ繋げる．

ToAdjacencyMatrix[AddEdges[GridGraph[5, 5], Join[Table[ {5 i , 5 i - 4}, {i, 1, 5}] , Table[{i, 5 (5 - 1) + i}, {i, 1, 5}]]]]//MatrixForm

一般化しておこう．

$makeRBCGraph[n_Integer] := AddEdges[GridGraph[n, n], Join[Table[ {n i , n i - (n - 1)}, {i, 1, n}] , Table[{i, n (n - 1) + i}, {i, 1, n}]]]$

グラフができてしまえば，それに対して small world を作ることはたやすい．上で作成したモジュールを一部変更すればよい．

[Graphics:HTMLFiles/index_215.gif]

接続を書くときに周期境界条件を使わないので中央に入れ替えが集中した印象があるが，見た目の問題だけなので気にしないことにする．

次に，このグラフの統計的な性質について調べてみる．

と，たっぷりと入れ替えてしまうと孤立する頂点が出現してしまう．少しに押さえておくと，どうなるかというと，

直径はあまり変化がおきない．入れ替えがもう少し進まないと一番遠い部分がどうしても残ってしまうのかもしれない．そこで平均距離に注目してみることにする．幾つか確率を変えてみると，

と， p が大きくなると計算ができないものの，小さくなっていくようである．

Scale Free Network

次に Barabási and Albert が提案した Scale Free Network について，その作り方と統計的な性質を見てみることにしよう．Networkの形成は２つのステップからなる．
１．成長（Growth）：小さな数の頂点から始めて,時間ステップが1だけ過ぎるごとに m (≤)本の辺を持つ1個の頂点を追加してゆく．辺は既存の点と繋がるようにする．
２．選択結合（Preferential　Attachment）：その点の次数に比例した確率で繋がる．つまり，個々の点の次数をとするとi 番目の頂点に連結する確率は =/∑ である.
このアルゴリズムによって作られたネットワークが Scale Free Network と呼ばれるのは，次数 k を持つ点の存在確率 P(k) がべき乗則に従うことがわかっているためである．

上記のプログラムが初期ヴァージョンだったのだが，m< のときに，最初に作る頂点が孤立してしまい，その後孤立したままになってしまう．そこで，最初の個については完全グラフを作るように変更した．それが以下のヴァージョンである．

たとえば =3, m= $RowBox[{1, _, Cell[]}]$ を100世代まで作ってみると以下のようになる．

[Graphics:HTMLFiles/index_238.gif]

シーケンシャルに追加の様子を見る場合には以下のような感じ．

=m=3 を50世代の場合には以下のようになる．

[Graphics:HTMLFiles/index_253.gif]

プログラムの説明や解析結果については後日．簡単にメモをつけておくと，最終結果の平均経路は，

calcMeanPath[Last[reslst]]

7230/1751

%%//N

-Graphics -

で，計算できる．もしも，グラフの形成の過程でどのように値が変化してゆくのかを知りたいのならば，以下の手順で行えばよい．

pathlst = Map[calcMeanPath, reslst] ;

ListPlot[pathlst, PlotJoined→True, Frame→True] ;

[Graphics:HTMLFiles/index_261.gif]

100 ぐらいで飽和しているようなので，それを確かめておく．

グラフの表示

Scale Free Network を円周上の点を用いて表示したいというリクエストが来たので，ランダムに点を置くものと，等間隔に点を置くものを作ってみた．もうひとつパラメータを置いてオーバーロードさせようかとも思ったけれど，紹介するだけにとどめる．

ShowGraph[Last[scaleFreeNetwork2[50, 4, 3]]] ;

[Graphics:HTMLFiles/index_264.gif]

ShowGraph[Last[scaleFreeNetwork3[50, 3, 2]]] ;

[Graphics:HTMLFiles/index_267.gif]

Combinatorica でパーコレーション

うまく書けるかどうかわからないけど，できそうなことは書いてみよう．とりあえず， combinatorica を読み込むところは同じである．

<<DiscreteMath`Combinatorica`

おそらく使うのは Grid グラフぐらいなものだろう．作り方をもう一度確認しておくと，

(grph2d = GridGraph[10, 10])//ShowGraph

[Graphics:HTMLFiles/index_270.gif]

-Graphics -

(grph3d = GridGraph[10, 10, 10])//ShowGraph

[Graphics:HTMLFiles/index_273.gif]

-Graphics -

Conbinatorica 的に言えば，パーコレーションとは辺の生き残りゲームを行うことだ．

s = InduceSubgraph[grph2d, RandomKSubset[Range[100], 80]] ;

ShowGraph[s] ;

[Graphics:HTMLFiles/index_277.gif]

上のようにひとつの方法として，生き残る点の数を与えてランダムに点を選ぶことで，新しいグラフを作ることができる．生き残る点の数を変えると，徐々に連結されていく様子がわかる．いわゆるサイトパーコレーションである．下の図はダブルクリックすると徐々に接続する辺の数が増えてゆく．しかし， RandomKSubset を使っているので辺は固定ではない．

Table[ShowGraph[InduceSubgraph[grph2d, RandomKSubset[Range[100], i]]], {i, 10, 90, 5}] ;

繋がった辺を消すことなく増やしてゆくためには，点の数を追加してゆくような関数を作ればよいはずである．少し格好悪いけど，こんな関数を作ってみよう．

randomSwap[Range[20]]

{8, 13, 18, 12, 9, 6, 11, 15, 17, 16, 5, 14, 1, 4, 20, 2, 7, 3, 19, 10}

randomSwap[{a, b, c, d, e}]

{c, e, a, d, b}

(randomSwap[Range[200]]//Sort) == Range[200]

True

この関数を作ると以下のようにして，関数をでたらめに番号を増やしていくことができる．

randomSwap[Range[10]]//Table[#[[Range[i]]], {i, 10}] &

{{6}, {6, 5}, {6, 5, 1}, {6, 5, 1, 10}, {6, 5, 1, 10, 9}, {6, 5, 1, 10, 9, 2}, {6, 5, 1, 10, 9, 2, 8}, {6, 5, 1, 10, 9, 2, 8, 3}, {6, 5, 1, 10, 9, 2, 8, 3, 4}, {6, 5, 1, 10, 9, 2, 8, 3, 4, 7}}

下のグラフをクリックして，頂点が増えることで浸透が進んでゆく様子を確認できる．

sublst = randomSwap[Range[V[grph2d]]]//Table[#[[Range[i]]], {i, V[grph2d]}] & ;

glst = Map[InduceSubgraph[grph2d, #] &, sublst] ;

Map[ShowGraph, glst[[Range[1, V[grph2d], 5]]]] ;

上の図はクリックすればわかれば，１回の描画で頂点が５点ずつ増えていく．何をやっているのかは明白だろう．

同じようにしてボンドパーコレーションを考えられないだろうか？確率を与えればよいのならこんな方法でいいのかな．

(g = DeleteEdges[grph2d, RandomKSubset[Edges[grph2d], 54]])//ShowGraph ;

[Graphics:HTMLFiles/index_331.gif]

Length[Edges[grph2d]]

180

とりあえずパーコレーションと同じ状態のグラフを作ることはできそうである．では，これを解析するにはどうしたらよいのだろうか？サイトパーコレーションの場合には，明らかにConnectedQ では解析が十分ではない．なぜなら，ひとつでも孤立した点があると False を返すこの関数は，端から端まで繋がったグラフがあるのかを判断するパーコレーションの枠組みに合わないからだ．逆にボンドパーコレーションの場合には ConnectedQ でもいいような気がしていたり．微妙だ…

文献を漁ってみると，ゴールド＆トボニチク「計算物理学入門」にこんなことが書いてあった．「格子の大きさが有限だと，格子の端から端まで連結したことの定義には任意性がある．たとえば，連結した経路の定義として，(i)格子全体を横方向か縦方向のどちらかに連結した経路，(ii)ある定められた方向，例えば縦方向に連結した経路，あるいは，(iii)横方向と縦方向の両方に連結した経路，いずれを用いてもよい」だそうだ．じゃぁ，どれを使おうかという話になるんだが，結局システムサイズ無限大に外装してしまえば，どれも変わらないということのようである．

ということで ConnectedQ を用いてグラフの成長を調べてみると…

Map[ConnectedQ, glst]

なんとなくいい感じなんだけど，最初の True が… ConnectedQ の定義を修正しようかとも思ったのだが，いちいち修正してから使うのも大変なので最初が邪魔なときは Rest を使って最初だけ切り落としてみることにする．

((Cases[Map[ConnectedQ, glst], False]//Length) + 1)/(glst//Length)//N

0.78

という感じで，全体が繋がるまでの確率が計算できた．サイトパーコレーションの確率としては大きすぎ．幾つか試してみることにする．

Table[((Cases[Map[ConnectedQ, makeConnectionSequence[grph2d]], False]//Length) + 1)/(glst//Length)//N, {10}]

{0.78, 0.81, 0.84, 0.74, 0.94, 0.75, 0.61, 0.8, 0.81, 0.77}

境界の影響かもしれないということで，周期境界条件の場合を試してみる．

ShowGraph[GridGraph[10, 10], VertexNumber→True, PlotRange→All] ;

[Graphics:HTMLFiles/index_342.gif]

10 × 10 のシステムならば，横方向は 10 ( i -1) +1 と 10i を結び，縦方向には i と 10×(10-1)+i を結べば良さそうだ．

Join[Table[{10 (i - 1) + 1, 10 i}, {i, 10}], Table[{i, 10 (10 - 1) + i}, {i, 10}]]

{{1, 10}, {11, 20}, {21, 30}, {31, 40}, {41, 50}, {51, 60}, {61, 70}, {71, 80}, {81, 90}, {91, 100}, {1, 91}, {2, 92}, {3, 93}, {4, 94}, {5, 95}, {6, 96}, {7, 97}, {8, 98}, {9, 99}, {10, 100}}

g = AddEdges[GridGraph[10, 10], %23]

AddEdges[-Graph:<180, 100, Undirected>-, -Graphics -]

Map[ConnectedQ, makeConnectionSequence[g]]

makeConnectionSequence[ConnectedQ[AddEdges[-Graph:<180, 100, Undirected>-, -Graphics -]]]

全体の数に対する False の数の比が閾値になるんじゃないかと思われる．

Table[((Cases[Map[ConnectedQ, makeConnectionSequence[g]], False]//Length) + 1)/(glst//Length)//N, {20}]

{0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01}

少し小さくなったけど，まだ大きめ．計算時間がかかるのを覚悟してシステムサイズを大きくしてみると…

$makeRepeatedBoundaryGrid[n_Integer] := AddEdges[GridGraph[n, n], Join[Table[{n (i - 1) + 1, n i}, {i, n}], Table[{i, n (n - 1) + i}, {i, n}]]]$

$Table[((Cases[Map[ConnectedQ, makeConnectionSequence[makeRepeatedBoundaryGrid[20]]], False]//Length) + 1)/(20^2)//N, {10}]//Timing$

{142.797 Second, {0.7575, 0.7825, 0.84, 0.7525, 0.805, 0.6925, 0.7775, 0.9025, 0.7625, 0.8025}}

やはり大きめに出てしまう…もっとシステムサイズを大きくして外挿を試みるかとも思ったのだけど，仕事で使っている
コンピュータがずっと占有されてしまったので３日で中止した．後日，どこかで走らせてみることにする．

$(res = Table[Table[((Cases[Map[ConnectedQ, makeConnectionSequence[makeRepeatedBoundaryGrid[i]]], False]//Length) + 1)/(i^2)//N, {100}], {i, 5, 50, 5}]) ;//Timing$

$Aborted

さまざまな network 上での囚人のジレンマゲーム

「囚人のジレンマ」はゲーム理論ではもっとも重要なゲームである．この囚人のジレンマをある戦略に基づいて繰り返し実施したときに結果がどのように変わっていくのかを以下の要領でしらべてみることにする．戦略の選択の仕方などの詳細は Hauert and Doebeli を参照のこと．

プログラム作成のための簡単な覚書

Combinatorica のパッケージを読む．

<<DiscreteMath`Combinatorica`

いい加減なグラフを作ってそのグラフで動作確認してゆく．

rg2004 = RandomGraph[20, 0.4] ; ShowGraph[rg2004, VertexNumber→True] ;

[Graphics:HTMLFiles/index_358.gif]

adjacencies = ToAdjacencyLists[rg2004]

例えば５番目の vartex と連結していある相手を知りたい場合には５番目の要素を抽出すればよい．

adjacencies[[5]]

{10, 13, 14, 15, 16}

Thread を使ってふたつの連結している点を関数 f に送ろう．

Thread[f[5, adjacencies[[5]] ]]

{f[5, 10], f[5, 13], f[5, 14], f[5, 15], f[5, 16]}

各プレーヤーの現在の strategy を表すリストと payoff を計算する関数を作ると計算が可能になりそう．

はじめに strategy を表すリストを作る．0 ならば Cooperation ， 1 ならば Betray としよう．

strategy = Table[Random[Integer], {V[rg2004]}]

{1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1}

payoff を計算する関数は 0 と 1 の組み合わせが４つ必要なだけである．例えば cost-benefit ratio を r とすると，こんな感じになるだろう．

payoff[a_, b_, r_:1.] := If[a == 0, If[b == 0, 1, -r], If[b == 0, 1 + r, 0]]

Table[Table[payoff[i, j, r], {i, 0, 1}, {j, 0, 1}]//MatrixForm, {r, 0.1, 0.4, 0.1}]

{( {{1, -0.1}, {1.1, 0}} ), ( {{1, -0.2}, {1.2, 0}} ), ( {{1, -0.3}, {1.3, 0}} ), ( {{1, -0.4}, {1.4, 0}} )}

自分と相手の戦略から自分の利益を計算する関数を定義する．

calcBenefit[x_Integer, y_Integer] := payoff[strategy[[x]], strategy[[y]], 0.4]

連結するすべての相手とゲームをしたときの利益の計算例．

Apply[Plus, Thread[calcBenefit[5, adjacencies[[5]] ]] ]/Length[adjacencies[[5]] ]

0.56

すべての点についてそれを行ってみる．

benefitList = Table[Apply[Plus, Thread[calcBenefit[i, adjacencies[[i]] ]] ]/Length[adjacencies[[i]]], {i, V[rg2004]}]

{0.7, -0.2, -0.05, 0.4, 0.56, -0.4, -0.225, 0.525, 0.56, 0.525, 1.05, 0.7, 0.6, 0.466667, -0.12, -0.2, 0.933333, 0.16, 0.466667, 0.525}

となって，めでたく利益のリストを得ることに成功した．このリストを元にして戦略を変更することになる．

自分と連結している中からランダムに選んだ相手の利益を比べてみて，相手の方が利益が大きい場合，自分の戦略を変更することにしよう．自分の位置が与えられたときに，新しい戦略を返す関数を定義する．

degs = Degrees[rg2004]

{6, 7, 8, 7, 5, 2, 8, 8, 5, 8, 4, 8, 7, 9, 5, 7, 6, 10, 6, 8}

Table[Random[Integer, {1, degs[[5]]}], {50}]

{3, 3, 1, 1, 1, 4, 1, 1, 4, 2, 1, 2, 4, 4, 1, 5, 4, 4, 2, 3, 5, 5, 2, 1, 5, 4, 1, 3, 3, 4, 5, 4, 5, 2, 3, 5, 5, 3, 2, 2, 3, 4, 1, 4, 3, 5, 2, 5, 4, 4}