plot(1:max(degree(lcc)),degree.distribution(lcc)[-1],
     log='xy',col='#377EB8',pch=20,
     xlab='Degree', ylab='Number of Vertices', main='')


# plot and fit the power law distribution
# calculate degree
d = degree(graph, mode = "all")
dd = degree.distribution(graph, mode = "all", cumulative = FALSE)
degree = 1:max(d)
probability = dd[-1]
# delete blank values
nonzero.position = which(probability != 0)
probability = probability[nonzero.position]
degree = degree[nonzero.position]
reg = lm(log(probability) ~ log(degree))
cozf = coef(reg)
power.law.fit = function(x) exp(cozf[[1]] + cozf[[2]] * log(x))
alpha = -cozf[[2]]
print(paste("Alpha =", round(alpha, 3)))
# plot
print(d)
curve(power.law.fit, col = "red", add = T, n = length(d))


# fit_power_law(lcc)