Bayesian inference for power laws

Question

Let $G$ be a graph with $N$ nodes. Let $p(d_i)$ be the probability of node $i$ to have $d$ connections. If this follows a power-law:

$$ p(d_i) = \frac{d_i^\alpha}{\sum_{j=1}^{N} d_j^\alpha} $$

$\alpha$ is often estimated by Maximum Likelihood Estimation.

Is there a common Bayesian estimation of the posterior?

$$ p(\alpha | G) \propto \prod_{i=1}^N d_i^\alpha \times p(\alpha) $$

Answer 1

The MLE estimate $\hat{\alpha}$ could be used as location parameter of a Normal distribution with scale parameter $\sigma$ used as prior distribution of the scaling parameter. Then, such a prior distribution can be updated into a posterior via Metropolis-Hastings algorithm, i.e. a Markov Chain Monte Carlo method used to obtain a sequence of random samples from a probability distribution for which direct sampling is difficult.

On the left: posterior distribution of the scaling parameter; light blue dotted lines represent the 90% High Density Interval, whereas the red solid line indicates the MLE estimate, which coincides with the median of the posterior distribution itself. On the right: the MCMC trace of the Metropolis-Hastings algorithm with 50,000 iterations and a burn-in of 5,000.

Hereafter the R code:

bayesian.estimate <- function(x, q = 0.75, nsim = 50000, burnin = 5000) {
require(poweRlaw)
options(digits = 3)
x <- x[x > 0]
ox <- displ(x)
xmin <- quantile(x, q)
ox$setXmin(xmin)
alpha.mle.x <- estimate_pars(ox)$pars

likelihood <- function(distr, alpha, minim) {
  ll <- dpldis(distr, minim, alpha, log = T)
  return(sum(ll))
}

prior <- function(alpha, alpha.mle){
  alpha.prior <- dnorm(alpha, mean = alpha.mle, sd = 0.25, log = T)
  return(alpha.prior)
}

posterior <- function(distr, alpha, minim, alpha.mle){
  return(likelihood(distr, alpha, minim) + prior(alpha, alpha.mle))
}

proposal.function <- function(param) {
  return(rnorm(1, mean = param, sd = 0.1))
}

metropolis <- function(distr, minim, alpha.mle, startvalue, iterations) {
  chain <- array(dim = c(iterations + 1 , 1))
  chain[1,] <- startvalue
  for (i in 1:iterations) {
    proposal <- proposal.function(chain[i,])

    probab <- exp(posterior(distr, proposal, minim, alpha.mle) -
                  posterior(distr, chain[i,], minim, alpha.mle))
    if (runif(1) < probab) {
      chain[i + 1,] <- proposal
    } else {
      chain[i + 1,] <- chain[i,]
    }
  }
  return(chain)
}

# ALPHA X
startvalue <- c(alpha.mle.x)
chain <- metropolis(x, xmin, alpha.mle.x, startvalue, nsim)

burn.in <- burnin
acceptance.rate <- 1 - mean(duplicated(chain[-(1:burn.in),]))
acceptance.rate

par(mfrow = c(1,2))
hist(chain[-(1:burn.in),1], nclass = 100, border = 0, col = rgb(0,0,0,.15),
   main = expression(paste("Posterior of ", alpha[x])), xlab = "MLE estimate = red line")
abline(v = median(chain[-(1:burn.in),1]), col = "skyblue", lwd = 2)
abline(v = alpha.mle.x, col = "red", lwd = 2)
abline(v = quantile(chain[-(1:burn.in),1], 0.05), col = "skyblue", lwd = 2, lty = "dotted")
abline(v = quantile(chain[-(1:burn.in),1], 0.95), col = "skyblue", lwd = 2, lty = "dotted")

plot(chain[-(1:burn.in),1], type = "l",
   xlab = "MLE estimate = red line", ylab = expression(alpha[x]),
   main = "MCMC Trace")
abline(h = alpha.mle.x, col = "red", lwd = 2)

message("\nmedian alpha x = ", median(chain[-(1:burn.in),1]),
"\nupper alpha x = ", quantile(chain[-(1:burn.in),1], 0.95),
"\nlower alpha x = ", quantile(chain[-(1:burn.in),1], 0.05),
"\nHDI 90% alpha x = ", quantile(chain[-(1:burn.in),1], 0.95) - quantile(chain[-(1:burn.in),1], 0.05))

posterior.x <- chain[-(1:burn.in),1]