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 .

Loading