3
$\begingroup$

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

$\endgroup$
6
$\begingroup$

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.

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]
$\endgroup$
  • $\begingroup$ 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) $\endgroup$ – alberto Oct 23 '15 at 10:48
  • 1
    $\begingroup$ Unfortunately, I can't help with approximation methods, but maybe someone here is more prepared than me :) $\endgroup$ – stochazesthai Oct 23 '15 at 10:49
  • $\begingroup$ 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! :) $\endgroup$ – alberto Oct 23 '15 at 10:54
  • 1
    $\begingroup$ 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. $\endgroup$ – Xi'an Oct 23 '15 at 13:46
  • $\begingroup$ 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. $\endgroup$ – alberto Oct 23 '15 at 14:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.