復雜網絡與網絡博弈
本文省略大量預備知識，並以極簡化的方式介紹網絡博弈和演化動力學。
常見三種復雜網絡
復雜網絡本質是一種高連接屬性的拓撲結構。在網絡上進行博弈過程，網絡結構本身對博弈過程會產生影響。
常用的三種網絡有：ER隨機網絡（Erdős Rényi Random）、BA無標度網絡（Barabási Albert Scale-free）、小世界網絡（Small-world）。
ER網絡和SW網絡屬於均質網絡，ER網絡中每個節點的度均相同，收益信息在網絡上的傳導同樣快，用於研究合作行為在“人人平等”的情況下演化的過程。而BA無標度網絡屬於異質網絡，其特點是度分布滿足冪律特征，hub節點所傳導的信息比一般節點更多。SW網絡充分體現“六度分隔”現象，常用於研究合作行為在小世界現象下如何演化。
通過特定參數下的模擬，這三種拓撲結構基本如下圖：
所有圖的規模N = 50，ER網絡的連邊概率p = 0.01；SW網絡中每個節點的鄰居數量nei = 3，重連接概率為p = 0.2；SF無標度網絡的平均度<k> = 6。
網絡結構上，視每個節點為一個博弈主體，而將每條邊視為玩家間的直接聯系，這樣的種群被稱為結構種群。
具有收益矩陣的常見網絡博弈
網絡博弈多種多樣，我們能在不同的博弈類型中得到合作行為演化的不同軌跡。常見的博弈框架根據收益計算的不同基本分為具有收益矩陣的、非收益矩陣表達的。
這里討論三種具有單邊收益矩陣的兩人博弈，和一種非矩陣表達的博弈，分別為：
- 囚徒困境（PDG）
- 雪堆博弈（SDG）
- 捐贈博弈（DG）
- 最后通牒（UG）
PDG、SDG、DG 的單邊收益矩陣，UG 的收益計算形式分別表示為：
上圖中的參數 b 和 c 通常用來控制困境強度，以便更好的研究不同程度的困境將如何影響合作行為的演化，通常將 c 設置為 1，控制 b 值從大到小逐步變化。
其中最后通牒博弈中的 (\(p_i\), \(q_i\)) 代表第 i 個 player 的策略，一般習慣在 0-1 上取值，當 i 與鄰居j進行博弈時，\(p_i\) 表示個體 i 最多能提供給鄰居的金額，\(q_i\) 表示個體 i 最少能接受鄰居提供的金額。
上述四種博弈均只討論兩玩家的情形。
當重復博弈進行一輪時，個體 i 通過網絡連邊中的直接連接關系（鄰接矩陣中有連邊的部分），可以分別計算出它與所有鄰居 j 博弈所得的收益，進而通過加總全部收益得到累計收益（Total Payoff）。
定義 M 為博弈矩陣，f 為博弈所帶來的收益，S 表示個體選擇的策略，這里個體間的策略我們抽象地概括為親社會性的 “C” 策略、自私性的 “D” 策略。另外記 G 為結構種群。
則S可以表示為：
初始狀態下，每個個體的策略隨機給定 C or D，進一步的，定義在網絡 G(V,E) 上全部與個體 i 的鄰居集合為 \(\Omega_i\) 。則收益計算可以表示如下：
具體在計算時實際上是個體 i 的策略向量與個體 j 的策略向量通過收益矩陣 M 進行向量乘法運算，可以表達為：
然而分別對每一個體i反復迭代運算無形中增加計算負擔，於是根據原理，可以利用鄰接矩陣對上式重新定義，首先定義鄰接矩陣 A：
再記 \(\Phi_{ij}\) 為全部 \(S_i^T M S_j\) 構成的矩陣，則直接將 \(\Phi_{ij}\) 與鄰接矩陣做 Hadamard 積（對應位置相乘），i 與 j 不相連的項全部乘零消失，留下的項均為有連邊的收益計算結果，收益矩陣表示為：
自然的，個體 i 的累計收益表示為 F 的第 i 行或列求和即可：
根據上述步驟進行一輪博弈，計算出全部個體的累計收益向量后需要更新每個個體的初始策略。常見的動力學更新規則有：
- 模仿更新
- 從眾更新
- 期望更新
- 自反更新
這里的只采用 Fermi 更新，一種簡單的模仿更新。該模仿規則能夠根據個體收益差計算出個體 i 的策略互換概率 w，即：
其中 \(\kappa\) 可以理解為決策時產生的擾動，一般取很小的值，例如0.01。
Fermi 更新在通過概率 w 決定是否要與對手互換策略時，如果個體 i 有多個鄰居之間相連，則在滿足 w 條件時，需要先選擇其中一個鄰居進行策略互換，這里無論策略好壞（所帶來的收益高低），每個鄰居被抽到的概率均等。
Fermi 在更新過程中也可以遵循兩種不同的動力學規則，即同步更新或者異步更新:
同步更新是指一輪博弈結束后，所有 N 個節點需要順序遍歷，也就是每個節點都需要進行且只進行一次更新。
異步更新指的是系統總共需要迭代 N 次（節點個數次），每次更新的對象從 N 個個體中等概率抽樣，也就是最極端的情況下可能有節點從未更新過，或者有節點重復更新 N 次。
事實上進行博弈的過程中你會發現，隨着博弈重復進行下去，在特定的更新機制下復雜系統中個體策略將不斷變化，直到一定步長時，C 策略和 D 策略各占據一定比例不再改變（非理想狀態下應該圍繞某個比例小幅度波動）。
動力學數據生成實現
生成的數據主要是累計收益向量 \(Y = \sum{F_{N \cdot}}\)，以及每輪博弈的矩陣\(\Phi_{ij}\)，即有多少輪模擬就會有多少列Y，有多少個節點就會有多少個 \(\Phi_{ij}\)
原創，僅供自己回憶使用，請勿轉載
For R
gen_dynamic<-function(N, L, aver_k=6, net_type="BA", game_type="PDG"){
# Calculate the payoff between i and j
# which need Payoff Matrix（for others except SUG）
payoff<-function(si, sj, M){
return(t(as.matrix(si)) %*% M %*% as.matrix(sj))
}
# which no need Payoff Matrix（for SUG）
payoff_sug<-function(s_i, s_j, S){
p_i<-S[,s_i][1]; q_j<-S[,s_j][2]; p_j<-S[,s_j][1]; q_i<-S[,s_i][2]
if ((p_i >= q_j) & (p_j >= q_i)) {p_sug<-p_j + 1 - p_i}
if ((p_i >= q_j) & (p_j < q_i)) {p_sug<-1 - p_i}
if ((p_i < q_j) & (p_j >= q_i)) {p_sug<-p_j}
if ((p_i < q_j) & (p_j < q_i)) {p_sug<-0}
return(p_sug)
}
# Generate different Game structure except SUG
gen_PDG<-function(m.R, m.S, m.T, m.P){
return(matrix(c(m.R, m.S, m.T, m.P), 2, 2, byrow=T)) # PDG Payoff Matrix
}
gen_SDG<-function(m.r){
return(matrix(c(1, 1-m.r, 1+m.r, 0), 2, 2, byrow=T)) # SDG Payoff Matrix
}
gen_DG<-function(m.b, m.c){
return(matrix(c(m.b-m.c, -m.c, m.b, 0), 2, 2, byrow=T)) # DG Payoff Matrix
}
################################### Simulation of Dynamics ############################################
# ER and SW networks are transformed into degree control
p_or_m<-function(N, k){
return(round(k / (N-1)))
}
# Generate different Complex Networks
if(net_type == "BA"){G<-barabasi.game(n=N, power=1, m=aver_k, directed = F)} # BA Scale-Free
if(net_type == "ER"){G<-erdos.renyi.game(n=N, p.or.m=p_or_m(N, aver_k), type="gnp", directed = F)} # ER Random
if(net_type == "SW"){G<-sample_smallworld(dim=1, size=N, nei=aver_k/2, p=.0001)} # SW Regular
# Init Payoff Matrix
if(game_type == "PDG"){M<-gen_PDG(1, 0, 1.2, 0)} # R、S、T、P參數隨意設置
if(game_type == "SDG"){M<-gen_SDG(1.2)} # r參數隨意設置
if(game_type == "DG"){M<-gen_DG(2, 1)} # b、c參數隨意設置
# Init Strategy
if(game_type == "SUG"){
# for SUG
S<-as.data.frame(matrix(NA, 2, N))
for (i in 1:N) {
S[,i]<-c(runif(1), runif(1))
}
}else{
# For others
S<-as.data.frame(matrix(NA, 2, N))
for (i in 1:N) {
if(runif(1) >= .5){
S[,i]<-c(0, 1)
}else{
S[,i]<-c(1, 0)
}
}
}
# Record Adj
real_adj<-as.matrix(get.adjacency(G))
storage.mode(real_adj)<-"integer" # 因為Adj全是0-1，所以強行由double型儲存為integer型
# init X, Y
Y<-data.frame(matrix(NA, L, N)) # 每一列為一個時間步的累計收益
X<-list() # 為包含N個矩陣的列表（即分塊對角矩陣分塊存儲）
for (i in 1:N) {
X[[i]]<-matrix(NA, L, N)
}
# simulation process
for (l in 1:L) {
F_temp<-matrix(NA, N, N)
if(game_type == "SUG"){
for (u in 1:N) {
for (v in 1:N) {
F_temp[u, v]<-payoff_sug(u, v, S)
}
}
}else{
for (u in 1:N) {
for (v in 1:N) {
F_temp[u, v]<-payoff(S[,u], S[,v], M)
}
}
}
# Record Y and X
F_cumulate<-apply((real_adj*F_temp), 1, sum)
# F_cumulate<-apply((real_adj*F_temp), 2, sum)
Y[l, ]<-F_cumulate # The row of Y
for (n in 1:N) {
X[[n]][l,]<-F_temp[n,]
}
# Fermi Updating (Synchronization Update)
if(game_type == "SUG"){
for (i in 1:N) {
update_index<-sample(as.vector(neighbors(G, i, "all")), 1)
w<-1/(1 + exp((F_cumulate[i] - F_cumulate[update_index])/.1))
if(runif(1) >= w){
S[, i]<-S[, update_index]
S[, i]<-S[, i] + .05
}
}
}else{
for (i in 1:N) {
update_index<-sample(as.vector(neighbors(G, i, "all")), 1)
w<-1/(1 + exp((F_cumulate[i] - F_cumulate[update_index])/.1))
if(runif(1) >= w){
S[, i]<-S[, update_index]
}
}
}
}
# X block diag treatment
# X_diag_block<-as.matrix(bdiag(X)) #使用此代碼將X保存為分塊對角矩陣需要用到開頭備注的 “bigmemory” 包
return(list(Y = Y, X = X, real_adj = real_adj))
}
# dyn_vars<-gen_dynamic(100, 15, 6, "BA", "SUG")
# dyn_vars[["Y"]]; dyn_vars[["X"]]; dyn_vars[["real_adj"]]
For Python
import numpy as np
import math, random
import networkx as nx
# Gen Network
class Gen_Net:
def basf(n: int, k: int):
G_sf = nx.barabasi_albert_graph(n, k)
return G_sf
def ER(n: int, p: float):
G_er = nx.random_graphs.erdos_renyi_graph(n, p)
return G_er
def SW(n: int, k: int, p: float):
G_sw = nx.random_graphs.watts_strogatz_graph(n, k, p)
return G_sw
# Define Useful Tools
class Use_Tools:
def chunk(lst, size): # make selected array sort with rule size
return list(map(lambda x: lst[x * size: x * size + size], list(range(0, math.ceil(len(lst) / size)))))
def payoff(si, M, sj) -> float: # si*M*sj' and input is np.array
phi = np.matmul(np.matmul(si, M), sj.T)
return phi
# Define Game Structure
class Game:
def PDG(b: float):
M = np.array([[1, 0], [b, 0]]) # b = 1.2; M[0][0]=1; M[0][1]=0; M[1][0]=1.2; M[1][1]=0
return M
def SDG(r: float):
M = np.array([[1, 1 - r], [1 + r, 0]])
return M
# Initalize Strategies Set
G = Gen_Net.basf(100, 6) # average degree <k> = 6: k = 3
node_list = np.array(list(G.nodes)) + 1 # all nodes
edge_list = np.array(G.edges) + 1 # linkage's edges
nebmat = nx.to_numpy_matrix(G) # Adjacent Matrix of G
# np.savetxt('Adjacent_Matrix.csv', nebmat, delimiter = ',')
s_set = [] # every node which choose C or D with equal prob
for _ in range(len(node_list)):
if np.random.rand(1) > 1 / 2: s_set.append((0, 1))
else: s_set.append((1, 0)) # strategies set of #N nodes
# Monte-Carlo Simulation
L = 300; Phi = []
M = Game.PDG(1.2) # MC Rounds
Y = np.zeros([L, G.number_of_nodes()])
for _ in range(G.number_of_nodes()): # Store Phi
Phi.append(np.zeros([L, G.number_of_nodes()]))
for l in range(L):
for f in range(G.number_of_nodes()): # asynchronous updating
temp_payoff, neb_temp_payoff = [], []
focal_player = random.sample(list(node_list), 1)[0]
# print("\n", focal_player)
# find focal's neb and calculate payoff of focal player
rdn_neb = []
for i in range(len(edge_list)):
if edge_list[i][0] != focal_player:
if edge_list[i][1] == focal_player:
rdn_neb.append(edge_list[i][0])
temp_payoff.append(Use_Tools.payoff(np.array(s_set[edge_list[i][1] - 1]), M, np.array(s_set[edge_list[i][0] - 1])))
else:
rdn_neb.append(edge_list[i][1])
temp_payoff.append(Use_Tools.payoff(np.array(s_set[edge_list[i][0] - 1]), M, np.array(s_set[edge_list[i][1] - 1])))
# print("\n", rdn_neb)
# print("\n", temp_payoff)
f_i = np.sum(temp_payoff)
Y[l][focal_player - 1] = f_i
# find neb's neb and calculate payoff of player's neb
select_neb = random.sample(rdn_neb, 1)[0]
for i in range(len(edge_list)):
if edge_list[i][0] != select_neb:
if edge_list[i][1] == select_neb:
neb_temp_payoff.append(Use_Tools.payoff(np.array(s_set[edge_list[i][1] - 1]), M, np.array(s_set[edge_list[i][0] - 1])))
else:
neb_temp_payoff.append(Use_Tools.payoff(np.array(s_set[edge_list[i][0] - 1]), M, np.array(s_set[edge_list[i][1] - 1])))
# print("\n", select_neb)
# print("\n", neb_temp_payoff)
f_j = np.sum(neb_temp_payoff)
# Fermi updating
# fi_lambda = np.sum(nebmat[focal_player - 1]) # normlization
# fj_lambda = np.sum(nebmat[select_neb - 1])
# w = 1 / (1 + math.exp(((f_i / fi_lambda) - (f_j / fj_lambda)) / 0.1)) # k = .1
w = 1 / (1 + math.exp((f_i - f_j) / 0.1)) # k = .1
if np.random.rand(1) <= w: s_set[focal_player - 1] = s_set[select_neb - 1]
print("\n th {} Rounds of Fermi updating:\n {}".format(f, s_set))
# find Phi
for m in range(G.number_of_nodes()): # #(m) matrix
for c in range(G.number_of_nodes()): # #(col) of #(m)
Phi[m][l][c] = Use_Tools.payoff(np.array(s_set[m]), M, np.array(s_set[c]))
print("\n\n th {} Rounds of MC Simulation".format(l))
print("\n\n Y = \n {} \n\n Dim(Y) = {}".format(Y, Y.shape))
# check weather Y and Phi are correct
Y1_hat = np.matmul(np.matrix(Phi[0]),nebmat[:,0])
Y1 = Y[:,0].reshape(L, 1)
Y1 == Y1_hat