# Maximum likelihood estimator of the parameter of randomness in Watts and Strogatz's model (1998)

According to the paper

Menezes, M. B., Kim, S., & Huang, R. (2017). Constructing a Watts-Strogatz network from a small-world network with symmetric degree distribution. PloS one, 12(6), e0179120,

I wrote the R-code I attach to estimate via maximum likelihood the parameter of randomness $$p$$, based on Watts and Strogatz's model (1998), of a social network represented by its adjacency matrix. For more details about Watts and Strogatz's model, you can read

Watts, D. J., & Strogatz, S. H. (1998). Collective dynamics of ‘small-world’ networks. Nature, 393(6684), 440.

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

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 right 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){
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.
}


Now, I am going to build the sampling distribution of the maximum likelihood estimator (MLE) of the parameter of randomness $$p$$ by:

i) setting a value of $$p$$, e.g. $$p = 0.2$$,

ii) generating $$1000$$ social networks according to Watts and Strogatz's model, conditionally on a network density $$d$$, e.g. $$d = 0.1$$, a number $$N$$ of nodes, e.g. $$N = 100$$, and $$p = 0.2$$,

iii) estimating $$p$$ via maximum likelihood for each network generated, and

iv) obtaining the frequency distribution of these $$1000$$ maximum likelihood estimations of $$p$$.

p <- 0.2
d <- 0.1
N <- 100
MLE_p <- NULL
for(i in 1:1000){
K <- round(d * (N - 1)/2, 0)
size = N, nei = K, p = p)))
OptProblem <- nlminb(start = runif(1),

Following the histogram, the sampling distribution of the MLE of $$p$$ is not centered on $$p$$. According to statistical inference, this should not occur.