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

Question

我正在尝试实现一种算法来生成遵循Barabási-Albert (BA) model 的图形。在该模型下，度数分布服从幂律：

P(k) ~ k^-λ

指数 λ 应该等于 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
}

该图实际上似乎是无标度的，gamma=0.72±0.21 为 100 阶，gamma=0.68±0.24 为 10,000 阶，类似的结果各不相同参数m。但指数明显不同于预期的 gamma=3。

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

如果有人对 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);
  }
}

Answer 1

两件事可能会有所帮助：

sample_pa 获取指数的参数alpha = 3

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

幂律很难估计

仅在度数分布上运行 OLS/LM 会使您的指数接近 0 而不是 3（换句话说，被低估了）。相反，如果您使用具有高 xmin 的 igraph::power_law_fit 命令，您将得到接近 3 的答案。请查看 Aaron Clauset's page and publications 以获取有关估计幂律的更多信息。确实，您需要为每个度数分布估计一个最佳 x-min。

这里有一些代码会更好一些：

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)。改变它会改变你的答案。