我正在尝试实现一个在Barabási-Albert (BA)模型之后生成图的算法。在此模型下,度分布遵循幂律:
P(k) ~ k^-λ
其中指数λ应该等于3。
为了简单起见,我将重点放在R代码上,在这里我使用igraph
函数。然而,我使用的是示例问题1 != 3,这似乎是一个广泛讨论的话题( 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
}
图看起来实际上是无标度的,给出了10,000阶的gamma=0.72±0.21和10,000阶的gamma=0.68±0.24,并且类似的结果改变了参数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);
}
}
发布于 2018-08-18 23:37:47
有两件事可能会有帮助:
获取指数sample_pa
参数的alpha = 3
实际上是power = 1
和m = 1
(查看维基百科文章中的定义与igraph::sample_pa文档-- power
论点并不意味着权力分配的程度)。
幂律很难估计。
只要在度分布上运行OLS/LM,指数就接近于0(换句话说,低估了)。相反,如果您使用igraph::power_law_fit
命令和高xmin
,您将得到接近3的答案。有关估计幂律的更多信息,请查看艾伦·克劳塞特的网页和出版物。实际上,您需要估计每个度分布的最佳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
)。改变它会改变你的答案。
https://stackoverflow.com/questions/51913069
复制相似问题