使用Barabási-Albert模型计算和理解无标度网络

Posted 2021-05-05

我正在尝试实现一个在Barabási-Albert (BA) model之后生成图的算法。在这种模式下，学位分布遵循幂律：

P（k）~k ^ -l

指数λ应该等于3。

为简单起见，我将重点关注R代码，我正在使用igraph函数。然而，我得到了λ！= 3的网络。似乎这是一个广泛涵盖的主题（example question 1，eq2，eq3），但我找不到令人满意的解决方案。

在R中，我使用igraph:::sample_pa函数生成遵循BA模型的图形。在下面的可重现的例子中，我设置了

# Initialize
set.seed(1234)
order = 100
v_degrees = vector()

for (i in 1:10000) {
  g <- sample_pa(order, power=3, m=8)

  # Get degree distribution
  d = degree(g, mode="all")
  dd = degree_distribution(g, mode="all", cumulative=FALSE)

  d = 1:max(d)
  probability = dd[-1]
  nonzero.position = which(probability !=0)
  probability = probability[nonzero.position]
  d = d[nonzero.position]

  # Fit power law distribution and get gamma exponent
  reg = lm (log(probability) ~ log(d))
  cozf = coef(reg)
  power.law.fit = function(x) exp(cozf[[1]] + cozf[[2]] * log(x))
  gamma = -cozf[[2]]
  v_degrees[i] = gamma
}

该图表实际上似乎没有标度，给出订单100的γ= 0.72±0.21和10,000的γ= 0.68±0.24，并且类似的结果改变了参数m。但是指数明显不同于预期的gamma = 3。

事实上，我试图用不同的语言实现这个模型（C ++，请看下面的代码），但是我得到的结果与低于3的指数相似。所以我想知道这是否是对BA模型的常见误解或者有什么问题。先前的计算符合幂律分布，与通常预期的相反，这是BA模型的正常行为。

如果有人对C ++感兴趣或更熟悉C ++，请参阅下面的附录。

附录：C ++代码为了理解下面的代码，假设一个对象类Graph和一个connect函数，该函数在作为参数传递的两个顶点之间创建了一个边。下面我给出两个相关函数BA_step和build_BA的代码。

BA_step

void Graph::BA_step (int ID, int m, std::vector<double>& freqs) {
  std::vector<int> connect_history;
  vertices.push_back(ID);

  // Connect node ID to a random node i with pi ~ ki / sum kj
  while (connect_history.size() < m) {
      double U (sample_prob()); // gets a value in the range [0,1)
      int index (freqs[freqs.size()-1]);
      for (int i(0); i<freqs.size(); ++i) {
          if (U<=freqs[i]/index && !is_in(connect_history, i)) { // is_in checks if i exists in connect_history
              connect(ID, i);
              connect_history.push_back(i);
              break;
          }
      }
  }

  // Update vector of absolute edge frequencies
  for (int i(0); i<connect_history.size(); ++i) {
      int index (connect_history[i]);
      for (int j(index); j<freqs.size(); ++j) {
          ++freqs[j];
      }
  }
  freqs.push_back(m+freqs[freqs.size()-1]);
  }

build_BA

void Graph::build_BA (int m0, int m) {

  // Initialization
  std::vector<double> cum_nedges;
  std::vector<int> connect_history;
  for (int ID(0); ID<m0; ++ID) {
      vertices.push_back(ID);
  }

  // Initial BA step
  vertices.push_back(m0);
  for (int i(0); i<m; ++i) {
      connect(m0, i);
      connect_history.push_back(i);
  }
  cum_nedges.push_back(1);
  for (int i(1); i<m; ++i) cum_nedges.push_back(cum_nedges[cum_nedges.size()-1]+1);
  cum_nedges.push_back(m+m);

  // BA model
  for (int ID(m0+1); ID<order; ++ID) {
      BA_step(ID, m, cum_nedges);
  }
}

答案

两件事可能有所帮助：

获得指数sample_pa的alpha = 3论据

真的是power = 1和m = 1（检查维基百科文章中的定义反对igraph :: sample_pa文档--- power论证并不意味着幂律分布的程度）。

幂律很难估计

只需在度分布上运行OLS / LM就可以得到一个接近于0的指数（换句话说，低估了）。相反，如果你使用高igraph::power_law_fit的xmin命令，你会得到更接近3的答案。检查Aaron Clauset's page and publications有关估算幂律的更多信息。真的，你需要估算每个度数分布的最佳x分钟。

这里有一些代码可以更好地工作：

library(igraph)
set.seed(1234)
order = 10000
v_degrees = vector()
for (i in 1:100) {
  g <- sample_pa(order, power = 1, m = 1)
  d <- degree(g, mode="all")
  v_degrees[i] <- fit_power_law(d, ceiling(mean(d))+100) %>% .$alpha
}
v_degrees %>% summary()
##   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  2.646   2.806   2.864   2.873   2.939   3.120

请注意，我组成了x-min使用（ceiling(mean(d))+100）。改变这将改变你的答案。