tests/test_bound.R

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

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

alpha <- 1
net <- simple_net(time_step = 10000, p_const = 0.1, alpha = alpha)
stats    <- get_statistics(as.PAFit_net(net), binning = FALSE,
                           only_PA = TRUE)
result_perfect <- empirical_estimate_perfect(net)
result <- empirical_estimate(stats$final_deg)

final_deg <- stats$final_deg
deg_vec   <- table(final_deg)
degree    <- as.integer(names(deg_vec))
non_zero  <- degree > 0
deg_vec   <- deg_vec[non_zero]
degree    <- degree[non_zero]
true      <- rep(1,length(degree))
for (jj in 1:length(degree)) {
  true[jj] <- sum(stats$n_tk[, degree[jj]] > 0 & stats$m_t > 0)
}

true <- true[-length(true)]
names(true) <- as.integer(degree[-length(degree)])

T        <- sum(final_deg) + length(final_deg) - 2
E        <- sum(final_deg)
N        <- length(final_deg)
p_hat    <- N/(E + N)
W_k      <- E - degree[-length(degree)]^2 + 1

W_k_positive           <- W_k
W_k_positive[W_k <= 0] <- min(W_k[W_k > 0])


names(W_k_positive) <- as.integer(degree[-length(degree)])


N_k_t_normalized       <- rep(0,length(degree) - 1)



for (i in 1:(length(deg_vec) - 1)) {
  N_k_t_normalized[i]  <- sum(deg_vec[(i + 1):length(deg_vec)])  
}
names(N_k_t_normalized) <- as.integer(degree[-length(degree)])
N_k_t_normalized    <- N_k_t_normalized / N_k_t_normalized[1] * W_k_positive[1]

#W_k_new <- pmax(W_k,N_k_t_normalized)
#W_k_new <- 1/2*(W_k_positive+N_k_t_normalized)
#W_k_new <- sqrt(W_k_positive*N_k_t_normalized)
W_k_new <- W_k; W_k_new[W_k <= 0] <- sqrt(W_k_positive[W_k <= 0] * N_k_t_normalized[W_k <= 0])

q <- 1 - p_hat
W_k_4 <- T - ((degree[-length(degree)] + 1)/2)^((1 + q)/q) + 1
W_k_5 <- T - ((degree[-length(degree)] + 1)/2)^((1 + q)/(2 * q)) + 1



names(W_k_new) <- as.integer(degree[-length(degree)])
names(W_k_4) <- as.integer(degree[-length(degree)])
names(W_k_5) <- as.integer(degree[-length(degree)])

y_min <- min(c(deg_vec,W_k_positive, N_k_t_normalized, true, W_k_new ))
y_max <- max(c(deg_vec,W_k_positive,N_k_t_normalized,true, W_k_new ))

plot(as.integer(names(true)),true, xlab = "Deg", ylab = "Appearance time",
     log = "xy", col = "black")
points(as.integer(names(W_k_positive)), W_k_positive, col = "red")
points(as.integer(names(N_k_t_normalized)), N_k_t_normalized, col = "green")
#points(as.integer(names(W_k_new)), W_k_new, col = "brown")
points(as.integer(names(W_k_4)), W_k_4, col = "purple")
points(as.integer(names(W_k_5)), W_k_5, col = "pink")