# Bayesian inference for power laws

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)$$

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]

• Cool! very nice, thanks! M-H was actually my first option. Do you know whether this is doable by analytic approximation methods? (Laplace, Variational Inference) – alberto Oct 23 '15 at 10:48
• Unfortunately, I can't help with approximation methods, but maybe someone here is more prepared than me :) – stochazesthai Oct 23 '15 at 10:49
• I'll leave the question open for a while, just in case (the think is I would like it to be fast since I want to update the estimation every time the graph evolves). Thank you very much again! :) – alberto Oct 23 '15 at 10:54
• Since the MLE is used in the "prior", this is an empirical Bayes resolution rather than a genuine Bayesian solution: the data is used twice. – Xi'an Oct 23 '15 at 13:46
• So now I understand what "empirical Bayes" is! :) @Xi'an it is emp. bayes unless your data is new, isn't it? In my case I'm thinking of learning the parameters for a set of past graphs and then "updating" them to make link predictions in new graphs while they evolve. – alberto Oct 23 '15 at 14:08