How can I implement Metropolis Hastings algorithm in an undirected connected graph in R?

I have posted a problem in Stack Mathematics regarding a Metropolis Hastings algorithm in graph which someone can read it here

(A code solution is written in the stack mathematics link but it is in the CoCalc and I do not know how to translate in R.)

In a nutshell the problem is: Consider a finite, undirected, connected graph 𝐺=(𝑉,𝐸) with vertex set 𝑉 and edge set 𝐸. We do not know the length of |𝑉| of the vertices of 𝐺, nor its structure. We can only see local information about 𝐺. E.g. if we are at a vertex 𝑥∈𝑉 we can see the neighbors of 𝑥, i.e. the vertices 𝑦∈𝑉 for which (𝑥,𝑦)∈𝐸, as well as how many neighbors 𝑥's neighbors have. Let us denote by 𝑑(𝑥) the degree of 𝑥∈𝑉 the number of neighbors of 𝑥.

Compute the transition probabilities of the chain {𝑋𝑛}𝑛∈𝑁 of the Metropolis-Hastings algorithm simulating the uniform distribution in 𝑉 using those of the random walk in 𝐺 as the proposal transition probabilities.

For simplicity let us assume that the graph has 5 vertices.

How can I write the MH algorithm in R for this specific problem or translate the Cocalc in the stack math answer ?

Solution

Edit, (i) ## CoCalc code in R (rather slow).

## CoCalc code in R.
library(igraph)

# g <- graph(c(1,2, 2,4, 4,1, 1,3), directed=FALSE) 
g <- graph_from_literal(A-B, B-D, D-A, A-C)

cur = 1
freq <- rep(0, vcount(g))
names(freq) <- as_ids(V(g))
nit <- 1E4

set.seed(1)
system.time({
for ( i in seq(nit) ) {
  neigh    <- V(g)[.nei(cur)]
  nextnode <- neigh[sample(length(neigh), 1)]
  if (runif(1) < degree(g, cur) / degree(g, nextnode) ){
    cur <- nextnode
  }
  freq[cur] <- freq[cur] + 1
}
}) 
freq <- freq / sum(freq)
freq

Output.

  user  system elapsed 
   4.63    0.05    5.13 

     A      B      D      C 
0.2471 0.2477 0.2463 0.2589

Regarding your point (ii) in Stack Mathematics. If in an (un)dirceted graph all cycles are even, then the graph is bipartite and aperiodic. As in this example. Degrees may be different.

library(igraph)
g <- graph_from_literal(A-X, X-B, B-Y, Y-A, A-Z)
degree(g)
mmm <- as.matrix(g[])
mmm <- mmm / rowSums(mmm)

## Calculate power serie.
mms <- mmm
mms <- mms %*% mmm; mms