tests/test_estimate_alpha.R

In thongphamthe/mcPAFit: Estimating Preferential Attachment from a Single Network Snapshot in Constant Time

rm(list = ls())
library("mcPAFit")

alpha <- 1
net <- simple_net(time_step = 10000, 
                  p_const = 0.01, alpha = alpha)
stats <- get_my_statistics(net$net)
result_perfect <- empirical_estimate_perfect(net$net)
result <- empirical_estimate(stats$final_deg)


final_deg <- stats$final_deg

result_square <- empirical_estimate_square(stats$final_deg)
sqrt(result$variance) / result$A

y_min <- min(result$A,result_square$A, result_perfect$A[as.character(result$k)], result_square$A_1, 
             result_square$A_2, result_square$A_3, result$lower_A)
y_max <- max(result$A,result_square$A, result_perfect$A[as.character(result$k)], result_square$A_1, 
             result_square$A_2, result_square$A_3, result$upper_A)
  
plot(result$k, result$A, log = "xy", pch = 20,
     #ylim = c(y_min,y_max),
     xlab = "Degree", ylab = "Attachment rate")
polygon(c(result$k,rev(result$k)) , 
        c(result$lower_A[as.character(result$k)],
          rev(result$upper_A[as.character(result$k)])), 
        col = rgb(1,0,0,0.1),border = NA)

lines(result$k,result$k^alpha)
points(result_square$k[as.character(result$k)], result_square$A[as.character(result$k)], pch = 20,col = "blue")
points(result_perfect$k[as.character(result$k)], 
       result_perfect$A[as.character(result$k)], pch = 20,col = "red")
points(result_square$k, result_square$A_1, pch = 20, col = "brown")


points(result_square$k, result_square$A_2, pch = 20, col = "yellow")
points(result_square$k, result_square$A_3, pch = 20, col = "green")
points(result_square$k, result_square$A_4, pch = 20, col = "purple")
points(result_square$k, result_square$A_5, pch = 20, col = "pink")

deg_vec <- table(stats$final_deg)
degree <- as.integer(names(deg_vec))
k_1 <- quantile(stats$final_deg,95/100)
abline(v = k_1)





result$alpha
result$ci
result_perfect$alpha
result_perfect$ci
#alpha_timeline



N_k_t <- result$N_k_t 
N_k_t <- N_k_t / N_k_t[1] * result_square$W_k[1]
y_min <- min(c(deg_vec,result_perfect$T_k, N_k_t))
y_max <- max(c(deg_vec,result_perfect$T_k,N_k_t))


plot(names(result_perfect$T_k),result_perfect$T_k, pch = 20,
     ylim = c(y_min, y_max), col = "black", log = "xy",
     xlab = "Degree", ylab = "Num. of appearances")
points(names(result_square$W_k), result_square$W_k, pch = 20, col = "red")

points(names(N_k_t),N_k_t, col = "blue")

points(as.integer(names(deg_vec[-1])), deg_vec[-1],
       col = "green")





#points(degree,deg_vec/PA, col = "green")