2018年9月14日金曜日

スモールワールドネス計算時にERグラフを何個作るか

 スモールワールドネス(以後SW)は次式で定義される。

  SW=(Cnet/Crand)/(Lnet/Lrand)
   Cnet: 扱っているネットワークのクラスター係数
   Crand: 扱っているネットワークと同じパラメタを持つERグラフのクラスター係数
   Lnet: 扱っているネットワークの平均最短距離長
   Lrand: 扱っているネットワークと同じパラメタを持つERグラフの平均最短経路長
   ERグラフ:Erdos-Renyi modelにより作られるネットワーク(ノード数nとリンク率p)

 Crand, Lrandが解析的に求められれば良いのだが不明。全通りの平均を求めれば確実だが、ノード数が多いと現実的な時間では終わらない。そのためErdos-Renyi modelでランダムネットワークを複数個作成してCとLを求め、その平均を用いることになる。
 では一体何個の平均をとればよいのか。
 まず1例をみてみる。個人的にこの先の未来で扱うデータはノード数10~100個程度ではないかと思うので、ひとまずノード数は80個としてみる。リンク率はとりあえずの0.5。比較するERグラフは100個でいいや。(適当)

-----------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
inf=np.inf

def shortest_path_length(lengthmat):
    n=lengthmat.shape[0]
    shortest_path_length=np.zeros((n, n))
    for i_start_node in range(n):
        i_length=np.zeros(n)+inf
        candidate=np.zeros(n)+inf
        candidate[i_start_node]=0
       
        while np.min(candidate)!=inf:
            i_length[np.argmin(candidate)]=np.min(candidate)
            candidate=i_length.reshape(n, 1)+lengthmat
            candidate=np.min(candidate, axis=0)
            candidate[i_length!=inf]=inf
       
        shortest_path_length[i_start_node]=i_length
       
    return shortest_path_length

def characteristic_path_length(lengthmat, glo=True):
    n=lengthmat.shape[0]
    shortest_path_length_mat=shortest_path_length(lengthmat)
    characteristic_path_length=np.sum(shortest_path_length_mat, axis=0)/(n-1)
   
    if glo==True:
        characteristic_path_length=np.mean(characteristic_path_length)
       
    return characteristic_path_length

def clustering_coefficient(weightedmat, glo=True):
    triangles=triangle(weightedmat)
    k=np.sum(wei2unwei(weightedmat), axis=0)
    k*=(k-1)
    k[k==0]=1
    clustering_coefficient=2*triangles/k
   
    if glo==True:
        return np.mean(clustering_coefficient)
   
    return clustering_coefficient

def wei2unwei(mat):
    n=mat.shape[0]
    adjacency=np.zeros((n, n))
    adjacency[mat>0]=1
    return adjacency   
   
def triangle(weightedmat):
    n=weightedmat.shape[0]
    triangles=np.zeros(n)
    for i in range(n):
        for j in range(n):
            for k in range(n):
                temp=weightedmat[i, j]*weightedmat[j, k]*weightedmat[k, i]
                triangles[i]+=pow(temp, 1/3)
 
    triangles/=2
    return triangles

def ERgraph(N, p):
    graph=np.zeros((N, N))
    i=0
    while i<N-1:
        j=i+1
        while j<N:
            if np.random.uniform(0, 1)<p:
                graph[i, j]=1
               
            j+=1
        i+=1
       
    graph=graph+graph.T
    return graph

n_nodes=80
p=0.5
cc_ergraph=[]
cpl_ergraph=[]

for i in range(100):
    temp=ERgraph(n_nodes, p)
    cc_ergraph.append(clustering_coefficient(temp))
    temp[temp==0]=inf
    cpl_ergraph.append(characteristic_path_length(temp))
   
plt.hist(cc_ergraph)
plt.xlabel('clustering coefficient')
plt.ylabel('frequency')
plt.title('cc_ergraph')
plt.show()
plt.hist(cpl_ergraph)
plt.xlabel('characteristic path length')
plt.ylabel('frequency')
plt.title('cpl_ergraph')
plt.show()
-----------------------------------------------


 ちょいと足りないくらいかな、という印象。
 次に、ノード数とER生成回数を変えてどうなるかを調べてみる。ノード数は50, 100個とし、それぞれに対してERグラフの生成回数を10, 30, 100, 200, 300, 500としたときの、Cと Lの値をヒストグラムで出してみる。また正規性の確認としてシャピロウィルク検定を行う。

-----------------------------------------------
N_NODES=[50, 100]
p=0.5
n_exams=np.array([10, 30, 100, 200, 300, 500])

for n_nodes in N_NODES:
    for n_exam in n_exams:
        cc_ergraph=[]
        cpl_ergraph=[]

        for i in range(n_exam):
            cpl_candidate=inf

            while cpl_candidate==inf:
                temp=ERgraph(n_nodes, p)
                cc_candidate=clustering_coefficient(temp)
                temp[temp==0]=inf
                cpl_candidate=characteristic_path_length(temp)

            cc_ergraph.append(cc_candidate)
            cpl_ergraph.append(cpl_candidate)

        fig=plt.figure(figsize=(16, 8))
        ax1=plt.subplot(121)
        ax2=plt.subplot(122)
       
        ax1.hist(cc_ergraph)
        ax1.set_xlabel('clustering coefficient')
        ax1.set_ylabel('frequency')
        ax1.set_title('cc_ergraph ('+str(n_nodes)+' nodes, '+str(n_exam)+' exams)')

        ax2.hist(cpl_ergraph)
        ax2.set_xlabel('characteristic path length')
        ax2.set_ylabel('frequency')
        ax2.set_title('cpl_ergraph ('+str(n_nodes)+' nodes, '+str(n_exam)+' exams)')
       
        plt.savefig(str(n_nodes)+'nodes_'+str(n_exam)+'exams.png')
       
        cc_temp=scipy.stats.shapiro(cc_ergraph)
        cpl_temp=scipy.stats.shapiro(cpl_ergraph)
        print('Shapiro-Wilk test (n_nodes='+str(n_nodes)+', n_exam='+str(n_exam)+')\ncc_ergraph: p='+str(cc_temp[1])+'\ncpl_ergraph: p='+str(cpl_temp[1])+'\n')
    
-----------------------------------------------
Shapiro-Wilk test (n_nodes=50, n_exam=10)
cc_ergraph: p=0.37558773159980774
cpl_ergraph: p=0.4100293517112732

Shapiro-Wilk test (n_nodes=50, n_exam=30)
cc_ergraph: p=0.14942386746406555
cpl_ergraph: p=0.8486078977584839

Shapiro-Wilk test (n_nodes=50, n_exam=100)
cc_ergraph: p=0.7921257615089417
cpl_ergraph: p=0.3289673328399658

Shapiro-Wilk test (n_nodes=50, n_exam=200)
cc_ergraph: p=0.039140429347753525
cpl_ergraph: p=0.4688810706138611

Shapiro-Wilk test (n_nodes=50, n_exam=300)
cc_ergraph: p=0.6537032723426819
cpl_ergraph: p=0.025576237589120865

Shapiro-Wilk test (n_nodes=50, n_exam=500)
cc_ergraph: p=0.22154365479946136
cpl_ergraph: p=0.8461586236953735

Shapiro-Wilk test (n_nodes=100, n_exam=10)
cc_ergraph: p=0.7797001600265503
cpl_ergraph: p=0.38541314005851746

Shapiro-Wilk test (n_nodes=100, n_exam=30)
cc_ergraph: p=0.6236131191253662
cpl_ergraph: p=0.6676943898200989

Shapiro-Wilk test (n_nodes=100, n_exam=100)
cc_ergraph: p=0.22231294214725494
cpl_ergraph: p=0.08100102841854095

Shapiro-Wilk test (n_nodes=100, n_exam=200)
cc_ergraph: p=0.09412982314825058
cpl_ergraph: p=0.5377574563026428

Shapiro-Wilk test (n_nodes=100, n_exam=300)
cc_ergraph: p=0.418567955493927
cpl_ergraph: p=0.4255199432373047

Shapiro-Wilk test (n_nodes=100, n_exam=500)
cc_ergraph: p=0.4211519658565521
cpl_ergraph: p=0.8973879814147949

-----------------------------------------------















 それぞれでCrand、Lrandがいくらになったのかを求めておくのを忘れてしまったのだが、分布がどうであれCrand、Lrand自体はn_examが小さくてもある程度の狭い範囲に収まっているようにみえる。なので最初にざっくりデータの傾向をみる時点ではそこまで大きな回数のERグラフを生成する必要はないようにも思う。2通りしかやってないので雰囲気でしか話せないが、ノード数と同じくらいでいいかな。厳密にデータ解析しようとしたら、最低でもノード数の2倍~4,5倍ってとこですか。
 その辺も含めて結局のところ、どの程度のブレ幅までなら許容できるかをその都度考えなきゃダメってことですかね。大げさな話、解析したい元のネットワークたちのCnetが0.3だったり0.7だったりとぶれてたら、Crandが0.45だろうが0.55だろうが解析はできそうですし。逆に0.49と0.51の違いを扱いたい場合には、かなりの回数が必要そう。

 それと、ここに載せてない試行の中で、n_examが大きいにも関わらずシャピロウィルク検定で正規性がないとされたことがままあった。そういうときでもヒストグラム的にはとんがりコーンになっていたので、今回の目的に対してシャピロウィルク検定はいまいちよくないのかもしれない。てっきり回数稼げば正規分布になるだろうと思い込んでいたので、これについては近日中に考えようと思う。

 また、もしかしてCrand、Lrandはノード数の影響を受けないかも?とも思い始めた。や、ノード数2つしか試してないのでとやかく言えないのだが、雰囲気で話しているのでお許し願いたい。上ではリンク率0.5で通してたので、こちらもいろいろ変えてみる必要がありそう。

 今後の課題
・ノード数と必要なERグラフの数の関係性
・上のリンク率による影響
・リンク率を変えたときのCrand、Lrandの変化をみる
・weighted networkへの応用(もともとweighted networkを扱いたかったのだが、weightedなランダムネットワークをどのように作るかを考えた結果、まずunweightedな場合についてちゃんと知ってからweightedな解析について検証しよう、というのが今回のモチベである)

0 件のコメント:

コメントを投稿