メモ：複雑ネットワーク解析

p30
次数相関のピアソン相関係数がなぜこの式になるのかがよくわかりませんでしたので、後日改めて考えて。まぁ前に考えた時もここわからなくて飛ばしたんだけどね。。。
なお次ページについても同様。分子の変形に苦戦。

p36
「線形代数の知識より」線形代数をちゃんとやっていない私目には理解不能なので、線形代数をやったら戻ってきてください。

p43
「v(i)とv(j)がこの確率よりも高い確率で隣接するならば」って何？

p61
「ネコとサルについては、領野と領野の結合を部分的に調べた様々な論文がある。これらの結果を統合して作られた脳全体の領野間のネットワークがあり、公開されている。」
面白そう

p72, p87
拡張グラフ or p=0のWSモデルにおけるクラスター係数
全ノードでクラスター係数は一緒なので、のクラスター係数はネットワークのクラスター係数に等しい。
ノードの前後それぞれ<k>/2個先のノードまでリンクを張っており、片側で
_(<k>/2) C _2
の三角形。
またノードをまたぐ三角形の数は
1+2+...+(<k>/2-1)
まとめると
C=[_(<k>/2) C _2 *2+{1+2+...+(<k>/2-1)}] / _<k> C _2
  =(3<k>-6)/(4<k>-4)

こう書くと、クラスター係数のCとcombinationのCがごちゃごちゃするな。

p77
「ネットワーク全体が連結となるのp>=log(N)/Nのときであることが知られている」
確かめたい

夏休みが始まる

　22歳にもなって夏休みだからと浮かれていてはいけないので、やることをやっていきたい。去年みたいに「休み明けの試験(CBT)にさえ通れば何もかもが許される」だと逆に学習意欲が下がる人間なので、少なくとも去年の夏よりは勉強に集中する。

夏休みは今日で3日目なのだけれど、まぁ昨日まではレポートやってましたので特に進捗はなく。7月いっぱいほとんど進捗なかったね。

本日を含めてぴったり5週間。その中に東医体が4日間と家族で旅行が3日間あるので、その7日間に何かは期待しないほうがいいかと思う。そんなわけで1週間に1冊ペースで読み物を進めていこうかな。

・複雑ネットワーク-基礎から応用まで-
・統計学入門(基礎統計学1)
・自然科学の統計学(基礎統計学3)

この3冊は確定。複雑ネットワークの本は半分くらい流し読みしたんだけどもう一回ちゃんと読む。できれば1週間かけたくない。まぁ東医体までの1週間は練習もあって疲れるだろうし、そんくらいのペースでもいいか。あーはやく引退したい。
予定があるのも7日間だけではないので、あと1冊分の枠はバッファーとして空けておく。もし余るようならその時はその時。

なお旅行もろもろの7日間も何か読もうということで
・コネクトーム:脳の配線はどのように「わたし」をつくり出すのか
を注文しました。神経科学における教養として、こういう一般向けの本を読むのもいいのかもなと思いまして。

これにより東医体→帰省→家族旅行の16日間の行程の荷物は、統計学入門(基礎統計学1)とコネクトーム本の2冊に絞る。基礎統計学1に使う時間が9日間になってしまうのが難点だけれど、

一週間

pythonで書いたGUIをlinuxで動かせるようにしたいという話。ターミナルからpython3 program.pyとして起動させた場合にはうまくいくのだが、デスクトップのランチャーから呼び出した場合にどうしても~/bin/bashのプログラムをwhichで探すところで正常に動作しなくなる。なぜなのか、ということで原因探してるのだが解決していない。

igraphは記事を更新し続けているのでもうすぐですね。土日は触らないかもなので来週中には。

やれたこと少ないな～。

igraphに取り組み始めた

先日ちょろっと触って放置しているRだが、ここにきてigraphの使い方を学ぶ必要ができたので、igraphを導入としてRに取り組むことになった。
マニュアルをとりあえずさらおうとしたのだがやはり全く頭に入ってこないので、ここに記録を残しておかないといけなさそうという次第。

【R基本】
・<-で代入は入力めんどい
・四則演算は普通
・==、!=、<=
・|で論理和、&で論理積、ふたつ重ねるとベクトルの全要素に対する一つのプール演算
・rm():変数の削除
・特殊変数 ちなみに後続の関数にぶち込むと真偽が返ってくる
　　NA: 足すとNAになる　未定義的な？　is.na(x)
　　NULL: 足すとnumeric(0)になる　NAとの使い分けはいまいちわからない is.null(x)
　　Inf, -Inf: 無限大　無限大だとis.finite(x)でFalse
　　NaN: not a number　0割る0とかやると出てくる　is.nan(0/0)==True

【ベクトル】
・c(v1, v2, ,,, ,vn)
・a:b でaからbの数
・rep(a, b): aをb回リピートしたベクトル オプションでtimesつけるとそれをさらに何回、eachをつけるとそれぞれ何回ずつ繰り返す
・seq(a, b, c): aからbまでcおきに
・length(vector): ベクトルの要素数
・dim(vector): ベクトルの次元数
・sum(vector): 合計
・mean(vector): 平均
・sd(vecotor): 標準偏差
・cor(vector1, vector2): 相関係数

【factor】いまいち使いどころがわからず実感わかないので割愛

【行列】
・matrix(a, b, c): b*cの行列　中身はすべてa
・cbind(v1, v2, ,,, , vn): 列ベクトルの連結
・rbind(v1, v2, ,,, , vn): 行ベクトルの連結
・matrix[2,]とかやると2行目、matrix[,2]で2列目
・matrix[1:2, 3:4]で1,2行目の3,4列目
・t(matrix)で転置
・忘れてたけど、配列の先頭indexは1
・行列の要素積は*、内積は%*%とこれまた入力めんどくせえな

【リスト】後日

【データフレーム】
・例
dfr1<-data.frame( ID=1:4,
                           FirstName=c("John", "Jim", "Jane", "Jill"),
                           Female=c(F,F,T,T),
                           Age=c(22,33,44,55)
                          )
・各要素にアクセスするときはdfr1$FirstNameといったように$で結ぶ

【条件分岐等】
for(i in 1:10)
{
    pass
}

while(x<10)
{
    pass
}

if(x==0) y<-0 else y<-y/x


【視覚化】後日

【ネットワーク生成】
・make_empty_graph(n): ノードn個、リンク0のネット(というか群w)
・make_full_graph(n): ノードn個、全結合のネット
・make_star(n): 1つのノードと(n-1)個のノードがリンクをはっているネット
・make_tree(n, children=x, mode='directed'): childrenは生える枝の数なので、真ん中はdegree=x, 末端はdegree=1, 中間はdegree=x+1ということになる。
・make_ring(n): n個のノードで円を作りましょう
・sample_gnm(n=100, m=40): Erdos-Renyi モデルは元はノード数nとリンク生成率pからなるみたいだが、この関数はmでリンクの数を設定する模様
・sample_smallworld(dim=2, size=10, nei=1, p=0.1): 先日作ったことから個人的二おなじみのWSモデル。
・sample_pa(n=100, power=1, m=1, directed=F): Barabasi-Albertpreferential attachment model。スケールフリー奴。
・rewire(net, each_edge(prob=0.1)): netのリンクをある確率で張り替える

【ネットワーク関連】
・library(igraph)でまずはパッケージを読み込みましょう
・net<-graph_from_data_frame(d=links, vertices=nodes, directed=T)の形でネットを作成。2つ(ノードとリンク)のデータフレームを統合してネットワークを形成するのがgraph_from_data_frame()ってことですかね。
・なお本文と順番が前後するが、逆にnetからリンクやノードのデータフレームを取りだすときは、as_data_frame(net, what="edges")でリンク、as_data_frame(net, what="vertices")。
・E(net):リンクの全羅列
・V(net):ノードの全羅列
・E(net)$type:リンクのタイプの羅列
・V(net)&blabla:リンクのblablaの羅列
・ようはgraph_from_data_frame()で統合してできたnetから、リンクとノードの中身にアクセスする方法ですな。
・plot(net, edge.arrow.size=.4, vertex.label=NA)みたいにして可視化
・simplify(net, remove.multiple=F, remove.loops=T)でネットワークの簡素化。remove.multiple=Tで同ノード間の同じ方向の矢印は和をとって一つにしますということ。remove.loops=Tで自分自身へのリンクを削除
・ただしsimplify()を用いると異なるタイプの値でも構わず和をとられてしまう
・as_edgelist(net, names=T): netから隣接リストを取りだす
・as_adjacency_matrix(net, attr="weight"): netから隣接行列を取りだす
・bipartite.projection(net)で2部グラフを形成
・なお2部グラフのリンクは、隣接行列をその転置したものとかけ合わせることで求めることができる

【plotにあたって】
・細かい設定は多すぎるので、これは覚えるものではないと信じていきます
・plot()のレイアウトを設定できる。種類はたくさん。
layout_as_star, layout_components, layout_in_circle, layout_nicely, layout_on_grid, layout_on_sphere, layout_randomly, layout_with_dh, layout_with_drl, layout_with_fr, layout_with_gem, layout_with_graphopt, layout_with_kk, layout_with_lgl, layout_with_mds等(マニュアルに記載のあるものはこれですべて。他にあるかは知らない。)
・力学的モデルに基づいた描画としてFruchterman-Reingoldモデルlayout_with_fr()、Kamad-Kawaiモデルはlayout_with_kk(net)
・par(mfrow=c(2,2), mar=c(0,0,0,0))とかやると、描画ゾーンが4分割される
・dev.off()で分割エンド
・通常plot()ではx=[-1, 1], y=[-1, 1]でrescaleされて表示されるので、自分で大きさを設定するときには、事前に設定したレイアウトlと拡大(縮小)率xを用いてplot(net, rescale=False, layout=l*x)とする


【分析にあたって】
上のようにしたって何も見えてこないので、特徴を見出していかなければならないという話。

・リンクの削除：delete_edges(net, edges)
例えば全てのリンクの重みでは多すぎる場合に平均以下の重みのリンクを消す場合には

cut.off<-mean(links$weight)
net.sp<-delete_edges(net, E(net)[weight<cut.off])

とすれば消せる。他にタイプ別に表示させる際なんかは

net.type1<-net-E(net)[E(net)$type==1]
net.type2<-net-E(net)[E(net)$type==2]

とすることでリンクを削除/分割することできる。また分割前のネットで決めたレイアウトを与えると、分割後のネットがそのレイアウトに従ったまま表示できる。(もはや日本語では何を言っているか自分でもわからない)

l<-layout_with_fr(net) #net1, net2の親玉net
plot(net.1, layout=l)　#分割後のnet1のプロットにnetのレイアウトを与える
plot(net.2, layout=l)　#分割後のnet2のプロットにnetのレイアウトを与える

【GUIでいじりたい】
tkid<-tkplot(net)　でおｋ
なおレイアウト情報はtkplot.getcoords(tkid)のようにして得られる。得られたレイアウト情報(ここではlとする)はplot(net, layout=l)として指定すればプロットできる。

【グラフの見方はhairball plotのみではない！！】
マニュアルにheatmap（隣接行列の数字を色の濃さにしたマス目みたいなやつ）の表示の仕方が記載されているので参照


【two mode networks】
色変えたりして一つの図に乗っけよう！

【ネットワーク解析】各特徴量の解釈については割愛　Ruminov等参照
・edge_density(net, loops=F)    #リンクの数/(n*(n-1))(directed networkの場合)
・reciprocity(net)    #directed networkにおいて両方向性のリンク
　なお、例えばsample_gnm(n=5, m=10, directed=T)で1つ両方向性のリンクが出た場合には、reciprocity(net)は0.2になる(2/10)
・transitivity(net, type="global")
　全体のtransitivity。type="local"にすると各ノードまわりのtransitivity。(つまりベクトルが返ってくる)
・triad_census(net)
　triangleを形成する全16モチーフごとの数
・diameter(net, directed=F, weights=NA)
　"""最長の最短距離"""
　なおネットワークが分断されてる場合には、リンクしている範囲で最長の最短距離
・get_diameter(net, directed=F)
　最長の最短距離の経路　複数個あっても一つしか表示しないっぽい
・degree(net, mode="all")
　各ノードのdegree。modeにはall, in, out, totalの4種があるのだが、allとtotalの違いがいまいちわからない。
・plot(net, vertex.size=degree(net, mode="all")*3)なんてやるとおしゃれになる(degreeが大きいほどノードを大きく表示する)
・degree_distribution(net, cumulative=F, mode="all")
　netの全ノードでの各degreeごとの頻度。cumulative=Tだと累積になる(始点1.0からだんだん下がっていって0.0になる)
　なので、累積のグラフ(正式名称がわからない)が出したいときは

deg_dist<-degree_distribution(net, cumulative=T, mode="all")
plot(x=0:max(degree(net, mode="all")), y=1-deg_dist)

・centr_degree(net, mode="in", normalized=T)
3つの値を返す。graph-level centralityが何かわからなかったけどこちら(http://www.stat.washington.edu/people/pdhoff/courses/567/Notes/l6_centrality.pdf)に書いてありました。式的にはdegreeなどの中心性を表す指標のばらつきが大きいとgraph-level centralityが大きくなるので、"中心"を担うノードがはっきりしているほど大きくなるような値ということでいいのかしら。
$res - 各ノードのdegree
$centralization - sum(max(res)-res)/(theoretical_max)ってところ
$theoretical_max - netと同じノード数で理論上最大のsum(max(res)-res)。netがdirected networkであれば2*(n-1)^2、undirected networkであればn*(n-1)ということでいいのだと思う。

・closeness(net)
各ノードから他全ての(n-1)個のノードへの最短距離の和の逆数
・centr_clo(net)
closeness centralityを求める。
$res - closeness(net)*n
$centralization -
$theoretical_max -

・eigen_centrality(net)
・centr_eigen(net)
・betweenness(net)
・edge_betweenness(net)
・centr_betw(net)　　　　　　　この辺はいったん保留で

・hub_score(net)
・authority_score(net)
Kleinberg's hub centrality scoreというものみたいです。不勉強なのでまだ理解できないのでそれだけ書いておいて次にいきます。

・最短距離の類
mean_distance(net, directed=F): 全2ノード間の最短距離の平均、リンクの逆行あり
mean_distance(net, directed=T): 全2ノード間の最短距離の平均、リンクの逆行なし
distances(net): 全2ノード間の最短距離をn*nの行列で返す
distances(net, weights=NA): 重さ無視
いろいろオプションもある。あるノードから他のノードへの最短距離とかももちろん取り出せる。最短経路を色つけるとかも。そういう時は調べながら書こう、うん。具体例はmanualにいろいろ書いてあるので、必要になったときに自分のタスクでやろうね、うん。

・cocitation(net)
何を求めているかさっぱりわかりません。

・cliques(net)
ネットワークの中で全結合してるサブネットワークを探す関数
・largest_cliques(net)
ネットワークの中で最大のノード数の全結合サブネットワークを返す関数

・クラスタリングの話
いろいろなアルゴリズムがあるのでそれぞれについては今後調べるということでどうか一つ(?)

ceb<-cluster_edge_betweenness(net)
clb<-cluster_label_prop(net)
cfg<-cluster_fast_greedy(net)

出力されるのは「どのように分割するか」という情報であってネットワークではないので、例えば上のように計算した後にcluster_edge_betweennessに基づいてplotする際には以下のようにする。

ceb<-cluster_edge_betweenness(net)
plot(ceb, net)

その他便利グッズ
length(ceb)     #グループ数
membership(ceb)     #各ノードの所属グループ番号
modularity(ceb)
clossing(ceb, net)    #各エッジがinner-group(FALSE)かinter-group(TRUE)か

はて毎日更新するなど言っていたが

先週一週間はちょっと病院実習で病んでたので無理でしたね～
事実進捗はほぼゼロだったので更新する内容も無かったのだけれど。

先週の研究室では2つのフリーソフトの遊び方もとい使い方を学びました。

そんなことより、必要に迫られたので金曜日からTkinterに取り組んでおります。オブジェクト指向プログラミング自体初めてなのですが、最初に取り組むにはちょうどいい難易度くらいな程度のモノなので、この3連休はそれに取り組んでました。

や、昨日は一日卓球してたんだけども。1か月ぶりの真夏に2部練すると吐く、ということがわかりました。

明日から平日ですが、4日間なので幾分マシですね。気張っていきましょう。

無駄なんかないさ～

Watts-Strogatsモデルは結局以下の形で幕を下ろします。
---------------------------------------------------------
import numpy as np

def mk_WSmodel(N, k_half, p, matrix):
    list=reorder(mk_extended_cycle(N, k_half, matrix=False))
    change_index=np.random.choice(list.shape[0], np.int(list.shape[0]*p), replace=False)

    for index in change_index:
        match=match_list(list, index, N)
        list[index, 1]=match[np.random.choice(match.shape[0], 1)]
   
    if matrix:
        list=list2matrix(list)
   
    return list

def mk_extended_cycle(N, k_half, matrix):
    list=np.zeros(2).reshape(1, 2)
    for i in range(N):
        j=1
        while j<=k_half:
            list=np.r_[list, np.array([i, i+j]).reshape(1, 2)]
            j=j+1
    list=list[1:, :]%N

    if matrix:
        list=list2matrix(list)

    return list

def reorder(list):
    for i in range(list.shape[0]):
        if zero_one()==1:
            temp=list[i, 0]
            list[i, 0]=list[i, 1]
            list[i, 1]=temp
    return list

def match_list(list, index, N):
    match=np.ones(N)
    match[np.int(list[index, 0])]=0
    for i in range(list.shape[0]):
        if list[i, 0]==list[index, 0]:
            match[np.int(list[i, 1])]=0
        elif list[i, 1]==list[index, 0]:
            match[np.int(list[i, 0])]=0
    return np.arange(0, N, 1)[match==1]

def zero_one():
    return np.random.choice(2, 1)

def list2matrix(list, N):
    link=np.zeros(N*N).reshape(N, N)
    for i in range(list.shape[0]):
        link[np.int(list[i, 0]), np.int(list[i, 1])]=1
        link[np.int(list[i, 1]), np.int(list[i, 0])]=1
    return link

def cal_L(list, N):
    mean_length=np.zeros(N)
    for i in range(N):
        mean_length[i]=np.sum(dykstra_neighborlist(list, N, i))/(N-1)

    L=np.sum(mean_length)/N
    return L

def cal_C(list, N):
    triangle=np.zeros(N)
    link=list2matrix(list, N)
    for i in range(N):
        neighbor_list=extract_neighborlist(list, i)
        k=neighbor_list.shape[0]
        if k>1:
            for x in range(k):
                y=x+1
                while y<k:
                    if link[np.int(neighbor_list[np.int(x), 1]), np.int(neighbor_list[np.int(y), 1])]==1:
                        triangle[i]=triangle[i]+2.0
                    y=y+1
       
            triangle[i]=triangle[i]/(k*(k-1))

    C=np.float(np.sum(triangle))/N
    return C


def dykstra_neighborlist(list, N, start):
    end=np.zeros(N)+inf
    end[start]=0

    while np.sum(end)==inf:
        candidate=np.zeros(N)+inf
   
        for i in range(N):
            if end[i]!=inf:
                neighbor_list=extract_neighborlist(list, i)
                for j in range(neighbor_list.shape[0]):
                    if end[np.int(neighbor_list[j, 1])]==inf:
                        length=end[i]+1
                        if candidate[np.int(neighbor_list[j, 1])]>length:
                            candidate[np.int(neighbor_list[j,1])]=length
                       
        candidate[end!=inf]=inf
        end[np.argmin(candidate)]=np.min(candidate)

    return end

def extract_neighborlist(list, i):
    neighbor1=list[list[:, 0]==i]
    neighbor2=list[list[:, 1]==i]
    for j in range(neighbor2.shape[0]):
        temp=neighbor2[j, 0]
        neighbor2[j, 0]=neighbor2[j, 1]
        neighbor2[j, 1]=temp

    neighbor_list=np.r_[neighbor1, neighbor2]

    return neighbor_list

昨日の続きやってる場合ではない

昨日のWSモデルの検証ですが平均距離についてサテ室のパソコンでさえ無限に時間がかかりそうなので、Nを減らして検証することにする。とにかくスモールワールド性ができているかが確認できればいいのです。

ということN=100にしたところ、平均距離を計算するのにサテ室のPCだと18秒、自分のPCだと200秒くらいかかるのでやはりサテ室で作業するのがはかどりますな。これなら現実的な時間である。
次にクラスター係数を求めようということでクラスター係数を求める関数。
-----------------------------------------
def cal_C(list, N):
    triangle=np.zeros(N)
    link=list2matrix(list, N)
    for i in range(N):
        neighbor_list=extract_neighborlist(list, i)
        k=neighbor_list.shape[0]
        if k>1:
            for x in range(neighbor_list.shape[0]):
                y=x
                while y<neighbor_list.shape[0]:
                    if link[np.int(neighbor_list[np.int(x), 1]), np.int(neighbor_list[np.int(y), 1])]==1:
                        triangle[i]=triangle[i]+2.0
                    y=y+1
         
            triangle[i]=triangle[i]/(k*(k-1))
 
    C=np.float(np.sum(triangle))/N
    return C
------------------------------------------
これを1回計算するのに、サテPCで現在70分経過。。。
非常に厳しい。とりあえず待って統計赤本読んでゐる。


そして150分が経過。
もう嫌になったので'python complex network analysis'でググったところNetworkXというパッケージがヒット。早速pipでインストールしようとしたらもうすでに入ってたよ！！水曜日の報告会までにちょろちょろと触って感触だけ確かめておこう。


(7/10追記)修正しました。どうりで時間かかるわけだよ！

Watts-Strogatzモデルのネットワーク検証中

とりあえずスモールワールドになってるかなーと調べるために、参考書と同じ設定でネットワークを作成して、平均距離とクラスター係数が同じくらいになるかどうかを調べてみようということで、N=1000、<k>=6(k_half=3)、p=0.05で平均距離を求めようと下のような関数を書いて計算させて現在1時間経過。。。
参考書だと1000回の平均出してるんだけど、私のPCでは厳しそうである。しかもクラスター係数も絶対時間かかるよなぁ。。。。。。
脳灰白質のROIならN=100ちょっとくらいだけど、SNSの解析とかするならN=1000くらいならすぐ到達しそうだよなーと思うと厳しい世界だ。

明日朝一でサテ室のパソコンで回してみるとする。以前何だかの計算をサテPCにやらせてみたらめちゃくちゃ速かったので。そして回してる間に明日はBAモデルもやってみようか。

----------------------------------------------------------------
def cal_L(list, N):
    mean_length=np.zeros(N)
    for i in range(N):
        mean_length[i]=np.sum(dykstra_neighborlist(list, N, i))/(N-1)
 
    L=np.sum(mean_length)/N
    return L


def dykstra_neighborlist(list, N, start):
    end=np.zeros(N)+inf
    end[start]=0
 
    while np.sum(end)==inf:
        candidate=np.zeros(N)+inf
     
        for i in range(N):
            if end[i]!=inf:
                neighbor_list=extract_neighborlist(list, i)
                for j in range(neighbor_list.shape[0]):
                    if end[np.int(neighbor_list[j, 1])]==inf:
                        length=end[i]+1
                        if candidate[np.int(neighbor_list[j, 1])]>length:
                            candidate[np.int(neighbor_list[j,1])]=length
                         
        candidate[end!=inf]=inf
        end[np.argmin(candidate)]=np.min(candidate)
 
    return end

def extract_neighborlist(list, i):
    neighbor1=list[list[:, 0]==i]
    neighbor2=list[list[:, 1]==i]
    for j in range(neighbor2.shape[0]):
        temp=neighbor2[j, 0]
        neighbor2[j, 0]=neighbor2[j, 1]
        neighbor2[j, 1]=temp
 
    neighbor_list=np.r_[neighbor1, neighbor2]
 
    return neighbor_list

Watts-Strogatzモデルのネットワーク作成

土曜日あるある：平日に寝なさ過ぎて土曜日は昼まで寝ている。

【ネットワーク解析】
np.random.choice()：replace=Falseで重複なし

今日はWatts-Strogatzモデルのネットワークを作成する関数。
ちゃんと作れているのかの検証が出来ないのでなんともいえないのだが、とりあえず何か隣接行列っぽいものを出力してくれる。Pythonでネットワークを可視化してくれるやつ絶対あるだろうし、その手のパッケージを探したほうがいいんだろう。いいんだろうけども、そんなん絶対WSモデル作る奴なんてあるだろうし、そういう意味で調べたら負けだと思っているが無駄である。
--------------------------------------------
import numpy as np

def mk_WSmodel(N, k_half, p, matrix):
    list=reorder(mk_extended_cycle(N, k_half, matrix=False))
    change_index=np.random.choice(list.shape[0], np.int(list.shape[0]*p), replace=False)
 
    for index in change_index:
        match=match_list(list, index, N)
        list[index, 1]=match[np.random.choice(match.shape[0], 1)]
     
    if matrix:
        list=list2matrix(list)
     
    return list

def mk_extended_cycle(N, k_half, matrix):
    list=np.zeros(2).reshape(1, 2)
    for i in range(N):
        j=1
        while j<=k_half:
            list=np.r_[list, np.array([i, i+j]).reshape(1, 2)]
            j=j+1  
    list=list[1:, :]%N
 
    if matrix:
        list=list2matrix(list)

    return list

def reorder(list):
    for i in range(list.shape[0]):
        if zero_one()==1:
            temp=list[i, 0]
            list[i, 0]=list[i, 1]
            list[i, 1]=temp
    return list

def match_list(list, index, N):
    match=np.ones(N)
    match[np.int(list[index, 0])]=0
    for i in range(list.shape[0]):
        if list[i, 0]==list[index, 0]:
            match[np.int(list[i, 1])]=0
        elif list[i, 1]==list[index, 0]:
            match[np.int(list[i, 0])]=0
    return np.arange(0, N, 1)[match==1]

def zero_one():
    return np.random.choice(2, 1)

def list2matrix(list):
    link=np.zeros(N*N).reshape(N, N)
    for i in range(list.shape[0]):
        link[np.int(list[i, 0]), np.int(list[i, 1])]=1
        link[np.int(list[i, 1]), np.int(list[i, 0])]=1
    return link

2日坊主回避

【ネットワーク解析】
与えられたノード数Nと各ノード間にリンクが出来る確率pから作成されたランダムネットワークを返す関数mk_random_network(N, p)
与えられた隣接行列が示すネットワークが分断されている場合に、各分団(?)間にリンクを追加して連結した隣接行列を返す関数connect_network(link)

数日間のまとめ。
---------------------------------------------
import numpy as np
inf=np.float('inf')

def count_network(link):
    group=np.zeros(link.shape[0])
    i=0
    while np.min(group)==0:
        i=i+1
        group[dykstra_linkmatrix(link, np.argmin(group))!=inf]=i
    return i

def dykstra_linkmatrix(link, start):
    end=np.zeros(link.shape[0])+inf
    pre_end=end.copy()
    end[start]=0
    while np.sum(end==pre_end)!=end.shape[0]:
        pre_end=end.copy()
        length=link+end.reshape(link.shape[0], 1)
        length[:, end!=inf]=inf
        if np.min(length)!=inf:
            end[np.argmin(length)%link.shape[0]]=np.min(length)
    return end

def mk_random_network(N, p):
    link=np.zeros(N*N).reshape(N, N)+inf
 
    for i in range(N):
        link[i, i]=inf

    i=1
    while i<N:
        j=0
        while j<i:
            if np.random.uniform(0.0, 1.0, 1)<=p:
                link[i, j]=1
                link[j, i]=1
             
            j=j+1
     
        i=i+1
     
    return link

def connect_network(link):
    group=np.zeros(link.shape[0])
    i=0
    while np.min(group)==0:
        i=i+1
        group[dykstra_linkmatrix(link, np.argmin(group))!=inf]=i
 
    while np.sum(group)!=group.shape[0]:
        con_group=np.random.choice(np.int(np.max(group)), 2, replace=False)+1
        indexs_group1=np.arange(0, group.shape[0], 1)[group==con_group[0]]
        index1=indexs_group1[np.random.choice(indexs_group1.shape[0], 1)]
        indexs_group2=np.arange(0, group.shape[0], 1)[group==con_group[1]]
        index2=indexs_group2[np.random.choice(indexs_group2.shape[0], 1)]
     
        link[index1, index2]=1
        link[index2, index1]=1
     
        i=0
        group=np.zeros(link.shape[0])
        while np.min(group)==0:
            i=i+1
            group[dykstra_linkmatrix(link, np.argmin(group))!=inf]=i
     

    return link

1日坊主回避

今日は統計赤本と、統計のための行列代数上が届くので後で受け取りに行く。なお財布の中身が足りない模様。

【Rによる統計的検定と推定】
・パッケージのインスコ：utils:::menuInstallPkgs()

・箱ひげ図：boxplot(y ~ x, data)

・2つの平均の検定(t検定)：t.test(data, paired, alternative)
・3つ以上の平均の検定(一元配置分散分析)：oneway.test(y ~ x, data, var.equal)
・2因子以上の平均の検定(二元配置分散分析)：
・2つの分散の検定：var.test(data, alternative)
・3つ以上の分散の検定：bartlett.test(y ~ x, data)

【ネットワーク】
以前ダイクストラ法復習として記事書いたんですけど、ネットワークが分断されてるものに対応してなかったので、対応させた。
とりあえず隣接行列の方のみ。

import numpy as np
inf=np.float('inf')

def dykstra_linkmatrix(link, start):
 end=np.zeros(link.shape[0])+inf
 pre_end=end.copy()
 end[start]=0
 while np.sum(end==pre_end)!=end.shape[0]:
  pre_end=end.copy()
  length=link+end.reshape(link.shape[0], 1)
  length[:, end!=inf]=inf
  if np.min(length)!=inf:
   end[np.argmin(length)%link.shape[0]]=np.min(length)
 return end

早い時間に記事として書いて一日の中で更新していくスタイル

内科回ってる時、朝回診前後でテンプレ書いて、昼食時、夕食時の2回更新するようなことをしていたので、それにならい「進捗が生まれた時点で書き足していくスタイル」はなかなかいいのではないかと思った。ただしプロブレムリストみたいに書き連ねるとそれぞれ全部に進捗を生み出さなければいけない気がしてしまうのでそれはしない。

【Rによる統計的検定と推定】
はっきり言ってつまらないし、サンプルデータも無いのはどうなのか。
とりあえずRをインストールしてひたすらプログラム試しては次へ進んでいる。

・Rは変数の代入が矢印だということだけ知っていたのですが、イコールでも大丈夫なんすね。

・t検定：t.test(data, paired, alternative)

だいたい60ページ分くらい進んだ。すでにやる気がない。
実践Rのほうがやる気は起きそうだけど、あれ終わるには時間かかりすぎるんだよなぁ～。

やることが色々

　研究室に通いだしてから半年。この半年で色々学んだ気もしますがそれでもやっぱり全然足りなくて、研究室演習が始まるまでの今後半年で同じだけ能力値が上がったとしても、それでもやっぱり全然足りない。
　とはいえ色々サーベイをしたおかげで、短期的または長期的に必要になってくる知識が何なのかとか、やりたいけどやっていくのはキツそうなこととか、色々見えてきたのはこの半年での収穫ではないのかな。

　なお報告としては、Courseraはちゃんと11weeksやりとげまして修了。とはいえ5週目前後までの内容以外あまりはっきりと覚えてないですね。クロスバリデーションとかバリアント・バイアスの話などを知る点以外は、やはり浅くて結局自分でなんとかせなあかんのですよ。英語の勉強としても微妙でしたね。英語読むのは多少早くなったかも？（ほぼ英語字幕読んでただけ）

　これからの半年は大切にしないといけないとぼーっと考えていたわけなのですが、やはり毎日何を学んだのか記録していくことが必要だと思うわけです。例えば去ること1年半前になんとなく勉強していた複素関数の内容ではっきりと覚えてることって全くなくて、なんか特異点の周りをグルグルしてたなぁ～くらいで、これじゃ何も学んでなかったのと同じ。せめてもう少し具体的に何してたか思いだせるようにならないと、意味がない。
　なので明日からは毎日学んだことを記録していくようにする。ほんのちょっとだとしても。ただし書けないこともあるというのが難しいところである。書けないことの記録は残せないのである。ま、せいぜい数学とプログラミングどれくらいやったかと、どんな論文読んだかとかを書く程度になる気はするね。


今後の予定
・R
具体的に何を進めていくかというのはあまり決めていませんが、Rを学ぶついでに統計について学んでいくスタイルになるかと思う。まず「Rによる統計的検定と推定」(オーム社)あたりをさらっと終わらせた後、「実践 R 統計分析」(オーム社)あたりを進めるのかな。まぁググりながら使える用になればオッケーですよ。

・complex network
これについては現在「複雑ネットワーク：基礎から応用まで」(近代科学社)を読んでまして、これはRとかPythonあたりで色々なネットワーク解析ツール使う上で知っておけばいいかなー程度に考えている。理論的な部分よりツールを使う部分の方が優先な気はする。

・統計
赤本を今日ポチりました。Rで統計について触っていくのと並行して進めていく。統計は基礎数学シリーズのをやったんだけど、複素関数同様ほぼ忘れてしまっている。あの頃に学んだことって全部跡形もなく消えてしまっているので、そういうことが無いようにちゃんと勉強の痕跡を残していきたい。
なお「統計のための線形代数 上」(丸善出版)もポチった。これは研究室演習の期間にがっつり取り組んで短期間で終わらせないとダメな気がする。
今気づいたんだけど、一年半前に後回しにした線形代数と統計が、今になってピンポイントで必要になってくるの、さすがにひどくないすか？

・発表準備
抄録を提出したので、まぁまったり進める。

・英語
大学院入試にTOEFLが必要ということを考えるとあまり放っておけないよなぁと。TOEFLのための勉強というよりはその下地作りとして、英語を聞く機会と話す機会は作っていきたい。話す機会って難しいよなぁ。

・医学
ほんまむり。サテ室のビデオ受講同期見てると、俺は年末爆死するんじゃないかと思うわけですが、だとしても実習を蔑ろにして生み出した時間でQBを解くのはポリシーに反するので、真っ当に生きて行きたい。(なお隙間時間で研究室関連のことをやっていることについては棚上げ)

スケジュールとしてはまず今読んでいる複雑ネットワークの本を早いとこ読み終えて、夏休みも使いつつ統計の赤本を8月中に終わらせる。