R/empirical_estimate.R

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

empirical_estimate <- function(final_deg, outlier_prop = 95/100) {
  deg_vec  <- table(final_deg)
  degree   <- as.integer(names(deg_vec))
  non_zero <- degree > 0
  T        <- sum(final_deg) + length(final_deg) - 2
  E        <- sum(final_deg)
  deg_vec  <- deg_vec
  degree   <- degree
  W_k      <- E - degree + 1
  names(W_k)    <- degree 
  result        <- rep(0,length(deg_vec) - 1)
  names(result) <- as.character(as.integer(degree[1:(length(deg_vec) - 1)]))
  N_k_t         <- rep(0,length(result))
  variance      <- rep(0,length(result))
  
  for (i in 1:(length(deg_vec) - 1)) {
    N_k_t[i]  <- sum(deg_vec[(i + 1):length(deg_vec)])  
    result[i] <- sum(deg_vec[(i + 1):length(deg_vec)]) / (deg_vec[i]) / W_k[i]
    variance[i] <- result[i]^2 / sum(deg_vec[(i + 1):length(deg_vec)])
  }
  variance <- variance / result[1]^2
  result   <- result / result[1]
  sd_log   <- sqrt(variance/ result^2)
  lower_A  <- exp(log(result) - 2 * sd_log)
  upper_A  <- exp(log(result) + 2 * sd_log)
  
  degree        <- as.integer(names(result))  
  names(degree) <- as.character(as.integer(degree))
  names(N_k_t)  <- as.character(as.integer(degree))  
  
  names(variance) <- as.character(as.integer(degree))  
    
   while (TRUE) {
      thresh     <- quantile(final_deg, outlier_prop)
      if (length(degree[degree <= thresh]) <= 2) {
          if (outlier_prop < 1) {
              outlier_prop <- outlier_prop + 0.005  
          } else {
              stop("outlier_prop is greater than 1!");  
          }
      } else break;
  } 
  A          <- result[degree <= thresh]
  k          <- degree[degree <= thresh]
  non_zero   <- A > 0 & k + 1 > 0
  regress    <- lm(log(A[non_zero])~log(k[non_zero] + 1))
  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, outlier_prop = outlier_prop,
              thresh = thresh, variance = variance, lower_A = lower_A,
              upper_A = upper_A,
              alpha = alpha_empi, ci = ci, N_k_t = N_k_t,
              W_k = W_k))
}