R/empirical_estimate_perfect.R

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

empirical_estimate_perfect <- function(net, outlier_prop = 95/100, deg_thresh = sqrt(max(final_deg))) {
  stats    <- get_statistics(as.PAFit_net(net), binning = FALSE,
                             only_PA = TRUE)
  final_deg <- stats$final_deg
  deg_vec  <- table(final_deg)
  degree   <- as.integer(names(deg_vec))
  non_zero <- degree > 0
  T        <- sum(final_deg) + length(final_deg) - 2
  deg_vec  <- deg_vec[non_zero]
  degree   <- degree[non_zero]
  W_k      <- rep(1,length(degree))
  for (jj in 1:length(degree)) {
      W_k[jj] <- sum(stats$n_tk[, degree[jj]] > 0 & stats$m_t > 0)
  }
  result  <- rep(0,length(deg_vec) - 1)
  names(result) <- degree[1:(length(deg_vec) - 1)]
  
  for (i in 1:(length(deg_vec) - 1))
    result[i] <- sum(deg_vec[(i + 1):length(deg_vec)]) / deg_vec[i] / W_k[i]
  names(W_k) <- as.integer(degree)
  W_k <- W_k[1:(length(W_k) -1)]
  
  
  result        <- result / result[1]
  degree        <- as.integer(names(result))  
  names(degree) <- as.character(as.integer(degree)) 
  
  while (TRUE) {
    thresh     <- quantile(final_deg, outlier_prop)
    if (length(degree[degree <= thresh]) <= 2) {
      outlier_prop <- outlier_prop + 0.005  
    } else break;
  }
  
  #thresh     <- quantile(final_deg, outlier_prop)
  thresh     <- quantile(final_deg, 1)
  A          <- result[degree <= thresh]
  k          <- degree[degree <= thresh]
  non_zero   <- A > 0 & k > 0
  regress    <- lm(log(A[non_zero])~log(k[non_zero]))
  names(regress$coefficients) <- c("offset","Attachment exponent")
  
  alpha_empi <- regress$coefficients[2]
  
  
  res                    <- df.residual(regress)
  if (res > 0)
    ci         <- confint(regress,"Attachment exponent")
  else ci <- "N"
  
  return(list(A = result, k = degree, 
              alpha = alpha_empi, ci = ci,
              outlier_prop = outlier_prop,
              thresh = thresh,
              T_k = W_k))
}

# empirical_estimate <- function(deg_vec, G = 20) {
#   bin_result <- binning(deg_vec, G = G)
#   result     <- rep(0,G)
#   
#   
#   for (i in 1:(length(bin_result$bin_hist) - 1))
#     result[i] <- sum(bin_result$bin_hist[i:length(bin_result$bin_hist)]) / bin_result$bin_hist[i]
#   
#   return(list(result = result, bin = bin_result$bin, bin_hist =  bin_result$bin_hist, 
#               center_bin =  bin_result$center_bin, 
#               start =  bin_result$begin_deg, end =  bin_result$end_deg, G = G))
# }