1 #' @title Estimate rank probabilities with Markov Chains  2 #' @description Performs a probabilistic rank analysis based on an almost uniform  3 #' sample of possible rankings that preserve a partial ranking.  4 #' @param P P A partial ranking as matrix object calculated with [neighborhood_inclusion]  5 #' or [positional_dominance].  6 #' @param rp Integer indicating the number of samples to be drawn.  7 #' @details This function can be used instead of [exact_rank_prob]  8 #' if the number of elements in P is too large for an exact computation. As a rule of thumb,  9 #' the number of samples should be at least cubic in the number of elements in P.  10 #' See vignette("benchmarks",package="netrankr") for guidelines and benchmark results.  11 #' @return  12 #' \item{expected.rank}{Estimated expected ranks of nodes}  13 #' \item{relative.rank}{Matrix containing estimated relative rank probabilities:  14 #' \code{relative.rank[u,v]} is the probability that u is ranked lower than v.}  15 #'  16 #' @references Bubley, R. and Dyer, M., 1999. Faster random generation of linear extensions.  17 #' *Discrete Mathematics*, **201**(1):81-88  18 #'  19 #' @seealso [exact_rank_prob], [approx_rank_relative], [approx_rank_expected]  20 #' @author David Schoch  21 #' @examples  22 #' \dontrun{  23 #' data("florentine_m")  24 #' P <- neighborhood_inclusion(florentine_m)  25 #' res <- exact_rank_prob(P)  26 #' mcmc <- mcmc_rank_prob(P,rp = vcount(g)^3)  27 #'  28 #' # mean absolute error (expected ranks)  29 #' mean(abs(res$expected.rank-mcmc$expected.rank))  30 #' }  31 #' @export  32 mcmc_rank_prob <- function(P, rp = nrow(P)^3) {  33 2  n.full <- nrow(P)  34 2  MSE <- which((P + t(P)) == 2, arr.ind = T)  35 2  if (length(MSE) >= 1) {  36 0  MSE <- t(apply(MSE, 1, sort))  37 0  MSE <- MSE[!duplicated(MSE), ]  38 0  g <- igraph::graph.empty()  39 0  g <- igraph::add.vertices(g, nrow(P))  40 0  g <- igraph::add.edges(g, c(t(MSE)))  41 0  g <- igraph::as.undirected(g)  42 0  MSE <- igraph::clusters(g)$membership  43 0  equi <- which(duplicated(MSE))  44 0  P <- P[-equi, -equi]  45  } else {  46 2  MSE <- 1:nrow(P)  47  }  48 49 2  init.rank <- as.vector(igraph::topological.sort(igraph::graph_from_adjacency_matrix(P, "directed")))  50 2  res <- mcmc_rank(P, init.rank - 1, rp)  51 2  res$expected <- res$expected + 1  52 2  expected.full <- c(0, n.full)  53 2  rrp.full <- matrix(0, n.full, n.full)  54 2  for (i in sort(unique(MSE))) {  55 2  idx <- which(MSE == i)  56 2  if (length(idx) > 1) {  57 0  group.head <- i  58 0  rrp.full[idx, ] <- do.call(rbind, replicate(length(idx), res$rrp[group.head, MSE], simplify = FALSE))  59 2  } else if (length(idx) == 1) {  60 2  rrp.full[idx, ] <- res$rrp[i, MSE]  61  }  62  }  63 2  expected.full <- res$expected[MSE]  64 2  for (val in sort(unique(expected.full), decreasing = T)) {  65 2  idx <- which(expected.full == val)  66 2  expected.full[idx] <- expected.full[idx] + sum(duplicated(MSE[expected.full <= val]))  67  }  68 2  return(list(expected.rank = expected.full, relative.rank = rrp.full))  69 }  70

Read our documentation on viewing source code .