R/empirical_estimate_binning.R

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

empirical_estimate_binning <- function(final_deg, G = 10, deg_thresh = round(sqrt(max(final_deg)))) {
  deg.max     <- max(final_deg)
  deg_vec_old <- table(final_deg)
  degree_old  <- as.character(as.integer(names(deg_vec_old)))
  
  # binning all the degree after the degree threshold
  
  G           <- deg_thresh + 2
  bin_names   <- c(0:(deg_thresh),rep((deg_thresh + 1),deg.max - (deg_thresh + 1) + 1))
  center_bin  <- rep(0,G)
  times_bin   <- rep(0,G)
  n_k_bin     <- rep(0,G)
  m_k_bin     <- rep(0,G)
  deg_bin     <- rep(0,G)
  
  E <- sum(final_deg)
  
  
  for (i in 1:length(final_deg)) {
      if (final_deg[i] != deg.max) {

          bin_index             <- bin_names[final_deg[i] + 1] + 1
          
          n_k_bin[bin_index]    <- n_k_bin[bin_index] + 1
          deg_bin[bin_index]    <- deg_bin[bin_index] + final_deg[i]
          
          center_bin[bin_index] <- center_bin[bin_index] + log(final_deg[i])
          times_bin[bin_index]  <- times_bin[bin_index]  + 1
      }
  }
  

  for (i in 1:(length(degree_old)-1)) {
      bin_index             <- bin_names[as.integer(degree_old[i]) + 1] + 1
     
      m_k_bin[bin_index]    <- m_k_bin[bin_index] + sum(deg_vec_old[(i+1):length(deg_vec_old)])
  }
  
  
  
  non_zero   <- m_k_bin > 0 & n_k_bin > 0
  
  #center_bin <- exp(center_bin[non_zero] / times_bin[ non_zero])
  center_bin <-c(0:(deg_thresh + 1))[non_zero]
  #degree   <- deg_bin[non_zero] / n_k_bin[non_zero]
  

  degree    <- center_bin
  #W_k           <- (E - degree^2 + 1)
  #W_k[W_k <= 0] <- min(W_k[W_k >= 0])
  
  W_k        <- E - degree + 1
  result     <- m_k_bin[non_zero] / n_k_bin[non_zero] #/ W_k
  A          <- result
  k          <- degree
  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]
  
  result <- result / result[1]
  res                    <- df.residual(regress)
  if (res > 0) {
    ci         <- confint(regress,"Attachment exponent")
  }else ci <- "N"
  
  return(list(A = result, k = degree, m_k_bin = m_k_bin, 
              n_k_bin = n_k_bin, alpha = alpha_empi, ci = ci))
}