Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- install.packages("sna")
- install.packages("igraph")
- library(sna)
- library(igraph)
- # This is the Log-likelihood function in Meneses et al.'s article (2017).
- # This function receives as inputs:
- # p: parameter of randomness, based on Watts and Strogatz's model, and
- # AdjMat: adjacency matrix of the social network.
- loglike <- function(p, AdjMat){
- n <- nrow(AdjMat) # n: number of nodes
- Dens <- sna::gden(AdjMat) # Dens: network density of the social network.
- K <- round(Dens * (n - 1)/2, 0) # K: number of neighbors a node has to its rights side in the regular lattice before rewiring.
- # logProbs: vector of the logarithm of the probabilities of observing the degree centrality of each node in social network, based on Menezes et al.'s paper.
- logProbs <- log(apply(X = matrix(seq(n), ncol = 1),
- MARGIN = 1,
- FUN = function(j){
- m <- sum(AdjMat[, j])
- LowLim <- max(c(2 * K - m, 0))
- UppLim <- min(c(n - 1 - m, 2 * K))
- a <- 1/2 * p * (n - 1 - 2 * K)/n
- b <- p/n
- Factor1 <- dbinom(x = LowLim:UppLim, size = 2 * K, prob = a)
- Factor2 <- dbinom(x = m - 2 * K + LowLim:UppLim, size = (n - 2) * K, prob = b)
- sum(Factor1 * Factor2)
- }
- ))
- sum(logProbs) # output: loglikelihood of observed data.
- }
- # This part of the program attemps to build the sampling distribution of the maximum likelihood estimator (MLE) of the randomness parameter p, with 0 <= p <= 1.
- p <- 0.2 # True value of the randomness parameter.
- # This part of the program simulates 1000 times a social network of N nodes with density d and parameter of randomness p. For each social network, I obtain the MLE of p. Then, I build the sampling distribution of the MLE of p.
- d <- 0.1
- N <- 2000
- MLE_p <- NULL
- for(i in 1:100){
- K <- round(d * (N - 1)/2, 0)
- AdjMat <- as.matrix(igraph::get.adjacency(igraph::watts.strogatz.game(dim = 1,
- size = N, nei = K, p = p)))
- OptProblem <- nlminb(start = runif(1),
- objective = function(x) -loglike(p = x, AdjMat = AdjMat),
- lower = 0, upper = 1,
- control = list(trace = TRUE))
- MLE_p <- c(MLE_p, OptProblem$par)
- }
- hist(MLE_p, col = "gray60")
Add Comment
Please, Sign In to add comment