Generate degree-distributions from a beta-binomial distribution

Degree distributions can be generated for networks with N nodes by drawing from a shifted beta-binomial distribution, such that

\(d_i \sim \text{Beta-binomial}(N−2,\pi, \rho) + 1\)

where the mean degree is

\(k = \pi (N-1)\)

and the variance is

\(\sigma^2 = \pi (N-2)(1-\pi)(1+ (N-3) \rho)\)

where \(\rho = 0\) implies a standard binomial distribution with probability of success . The distribution is shifted up one value because a typical beta binomial distribution incudes some 0s, and the network generation algorithm, requires all nodes to have at least one connection.

The probability density function (pdf) defines the distribution. Below is the pdf for an exponential and the shifted beta-binomial, respectively

\(f(x) = P(X = x) = \lambda e^{\lambda x}\)

\(f(x) = P(X = x) = \binom n {x-1} \frac{\text{Beta}(x-1+\pi, n-(x-1))}{\text{Beta}(\pi, \rho)}\)

# dexp or dbetabinom defines the density function at a given value of x
par(mfrow = c(1, 2))
x <- seq(1, 15, 0.11)
plot(y = dexp(x, rate = 1), x = x, type = "l", ylab = "P(X = x)")
lines(y = dexp(x, rate = 0.5), x = x, lty = 2)
legend("topright", legend = c("rate = 1", "rate = 3"), 
    lty = c(1, 2), bty = "o")

x <- seq(1, 15,1) # discrete distribution
p1 <- dbetabinom(x-1, size = 15-2, prob = 0.5, rho = 0)
p2 <- dbetabinom(x-1, size = 15-2, prob = 0.5, rho = 0.1)
res <- rbind(p1, p2)
dimnames(res) <- list(c("rho = 0", "rho = 0.3"), x)
barplot(res, beside = TRUE, legend.text = TRUE, ylab = "P(X = x)")

Generate degree distributions

Networks can be generated using the algorithm of Viger and Latapy (2005), implemented in the degree.sequence.game function of igraph, which uniformly samples simple, connected graphs given the provided degree sequence.

In this example, I set the mean degree to 5. The +1 at the end shifts it. rbetabinom samples from the distribution.

# generate degree distribution
mean <- 5
num.populations <- 15
prob <- (mean - 1) / (num.populations - 2)
dseq <- rbetabinom(num.populations, num.populations - 2, prob = prob, rho = 0.1) + 1

g <- degree.sequence.game(dseq, method = "vl")
plot(g)

# edg <- get.edgelist(g)

How to estimate parameters from data

# Example data
data <- rgamma(20, rate = 1/10, shape = 2)

# function to calculate the negative log likelihood
betabinomNLL <- function(x, params){# params[1]=a, param[2]=b
    -sum(  ( log((choose(n = nval, k = x - 1))) + 
    lbeta(x - 1 + params[1], nval - (x - 1) + params[2]) - 
    lbeta(params[1], params[2]) )   )
}

#optim(fn = betabinomNLL, par = c(1, 1), x = round(data), 
#      method = "L-BFGS-B", lower = .01, upper = 500000)