R/query.R

In mcmcabn: Flexible Implementation of a Structural MCMC Sampler for DAGs

################################################################################ query.R --- Author : Gilles Kratzer Document created: 13/02/2019 -

##-------------------------------------------------------------------------
## Function to query mcmc samples generated by mcmcabn
##-------------------------------------------------------------------------

query <- function(mcmcabn = NULL, formula = NULL) {

    if (!(!is.matrix(formula) || !is.null(formula) || !inherits(formula,"formula"))) {
        stop("Formula statment should be given using either a matrix or a formula. Alternativelly, the argument could be null to return all individual arc support")
    }

    if (is.matrix(formula)) {

        n <- length(mcmcabn$scores)
        n.var <- length(mcmcabn$data.dist)
        m.array <- array(data = 0, dim = c(n.var, n.var, n))
        for (i in 1:n.var) {
            for (j in 1:n.var) {
                if (formula[i, j] == 1) {
                  m.array[i, j, ] <- 1
                }
                if (formula[i, j] == -1) {
                  m.array[i, j, ] <- -1
                }
            }
        }

        prob.1 <- which(x = m.array == 1, arr.ind = TRUE)
        prob.2 <- which(x = m.array == -1, arr.ind = TRUE)

        res.1 <- apply(matrix(data = mcmcabn$dags[prob.1], ncol = sum(m.1), byrow = TRUE), 1, prod)
        res.2 <- apply(matrix(data = mcmcabn$dags[prob.2], ncol = sum(m.2), byrow = TRUE), 1, prod)

        out <- sum(res.1 * res.2)/n
        return(out)
    }

    if (is.null(formula)) {
        out <- apply(mcmcabn$dags, 1:2, mean)
        colnames(out) <- rownames(out) <- names(mcmcabn$data.dist)
        return(out)

    }

    if (inherits(formula,"formula")) {

        f <- as.character(formula)

        if (!grepl("~", f[1], fixed = TRUE)) {
            stop("Formula specifications should start with a ~")
        }

        if (grepl("-", f[2], fixed = TRUE)) {

            formula <- sapply(strsplit(f[2], "[-]"), function(x) if (length(x) > 1)
                paste(x[1], paste(x[-1], collapse = "+"), sep = "-") else x)

            f <- strsplit(x = formula, split = "-", fixed = TRUE)
            f.1 <- paste0("~", f[[1]][1])
            f.2 <- paste0("~", f[[1]][2])

            m.1 <- formula.mcmcabn(f = as.formula(gsub(" ", "", unlist(f.1), fixed = TRUE)), name = names(mcmcabn$data.dist))
            m.2 <- formula.mcmcabn(f = as.formula(gsub(" ", "", unlist(f.2), fixed = TRUE)), name = names(mcmcabn$data.dist))

            n <- length(mcmcabn$scores)
            n.var <- length(mcmcabn$data.dist)
            m.array <- array(data = 0, dim = c(n.var, n.var, n))
            for (i in 1:n.var) {
                for (j in 1:n.var) {
                  if (m.1[i, j] == 1) {
                    m.array[i, j, ] <- 1
                  }
                  if (m.2[i, j] == 1) {
                    m.array[i, j, ] <- -1
                  }
                }
            }

            prob.1 <- which(x = m.array == 1, arr.ind = TRUE)
            prob.2 <- which(x = m.array == -1, arr.ind = TRUE)

            res.1 <- apply(matrix(data = mcmcabn$dags[prob.1], ncol = sum(m.1), byrow = TRUE), 1, prod)
            res.2 <- apply(matrix(data = mcmcabn$dags[prob.2], ncol = sum(m.2), byrow = TRUE), 1, prod)

            out <- sum(res.1 * res.2)/n
            return(out)

        } else {


            m <- formula.mcmcabn(f = formula, name = names(mcmcabn$data.dist))

            # creat array
            n <- length(mcmcabn$scores)
            n.var <- length(mcmcabn$data.dist)
            m.array <- array(data = 0, dim = c(n.var, n.var, n))
            for (i in 1:n.var) {
                for (j in 1:n.var) {
                  if (m[i, j] == 1) {
                    m.array[i, j, ] <- 1
                  }
                }
            }

            prob <- which(x = m.array == 1, arr.ind = TRUE)

            out <- sum(apply(matrix(data = mcmcabn$dags[prob], ncol = sum(m), byrow = TRUE), 1, prod))/n

            return(out)

        }
    }
}  #eof