inst/examples/helpfunctions_test.R

In kkholst/mcif: What the package does (short line)

Posdef <- function (n, ev = runif(n, 0, 10)){
    Z <- matrix(ncol=n, rnorm(n^2))
    decomp <- qr(Z)
    Q <- qr.Q(decomp)
    R <- qr.R(decomp)
    d <- diag(R)
    ph <- d / abs(d)
    O <- Q %*% diag(ph)
    Z <- t(O) %*% diag(ev) %*% O
    return(Z)
}

SigmaGen <- function(vcv, ncauses, old=FALSE){
    if (old==FALSE){
        if (ncauses==2){
            sigma88 <- matrix(c(1,0,0,0,0,0,0,0,
                                0,1,0,0,0,0,0,0,
                                0,0,1,0,0,0,0,0,
                                0,0,0,1,0,0,0,0,
                                0,0,0,0,vcv[1,1],vcv[2,1],vcv[3,1],vcv[4,1],
                                0,0,0,0,vcv[1,2],vcv[2,2],vcv[3,2],vcv[4,2],
                                0,0,0,0,vcv[1,3],vcv[2,3],vcv[3,3],vcv[4,3],
                                0,0,0,0,vcv[1,4],vcv[2,4],vcv[3,4],vcv[4,4]),
                              nrow=8,ncol=8)

            A <- matrix(c(1,0,0,0,0,0,
                          0,1,0,0,0,0,
                          0,0,1,0,0,0,
                          0,0,0,1,0,0,
                          1,0,1,0,0,0,
                          0,1,0,1,0,0,
                          0,0,0,0,1,0,
                          0,0,0,0,0,1),
                        nrow=6, ncol=8)
            sigma <- A%*%sigma88%*%t(A)
        }
        if (ncauses==3){
            sigma12 <- matrix(c(1,0,0,0,0,0,0,0,0,0,0,0,
                                0,1,0,0,0,0,0,0,0,0,0,0,
                                0,0,1,0,0,0,0,0,0,0,0,0,
                                0,0,0,1,0,0,0,0,0,0,0,0,
                                0,0,0,0,1,0,0,0,0,0,0,0,
                                0,0,0,0,0,1,0,0,0,0,0,0,
                                0,0,0,0,0,0,vcv[1,1],vcv[2,1],vcv[3,1],vcv[4,1],vcv[5,1],vcv[6,1],
                                0,0,0,0,0,0,vcv[1,2],vcv[2,2],vcv[3,2],vcv[4,2],vcv[5,2],vcv[6,2],
                                0,0,0,0,0,0,vcv[1,3],vcv[2,3],vcv[3,3],vcv[4,3],vcv[5,3],vcv[6,3],
                                0,0,0,0,0,0,vcv[1,4],vcv[2,4],vcv[3,4],vcv[4,4],vcv[5,4],vcv[6,4],
                                0,0,0,0,0,0,vcv[1,5],vcv[2,5],vcv[3,5],vcv[4,5],vcv[5,5],vcv[6,5],                                                 0,0,0,0,0,0,vcv[1,6],vcv[2,6],vcv[3,6],vcv[4,6],vcv[5,6],vcv[6,6]),
                              nrow=12,ncol=12)

            A <- matrix(c(1,0,0,0,0,0,0,0,0,
                          0,1,0,0,0,0,0,0,0,
                          0,0,1,0,0,0,0,0,0,
                          0,0,0,1,0,0,0,0,0,
                          0,0,0,0,1,0,0,0,0,
                          0,0,0,0,0,1,0,0,0,

                          1,0,0,1,0,0,0,0,0,
                          0,1,0,0,1,0,0,0,0,
                          0,0,1,0,0,1,0,0,0,
                          0,0,0,0,0,0,1,0,0,
                          0,0,0,0,0,0,0,1,0,
                          0,0,0,0,0,0,0,0,1),
                        nrow=9, ncol=12)
            sigma <- A%*%sigma12%*%t(A)
        }
        return(sigma)
    }
    if (old==TRUE){
        sigma88 <- matrix(c(1,0,0,0,0,0,0,0,
                            0,1,0,0,0,0,0,0,
                            0,0,1,0,0,0,0,0,
                            0,0,0,1,0,0,0,0,
                            0,0,0,0,vcv[1,1],vcv[2,1],vcv[3,1],vcv[4,1],
                            0,0,0,0,vcv[1,2],vcv[2,2],vcv[3,2],vcv[4,2],
                            0,0,0,0,vcv[1,3],vcv[2,3],vcv[3,3],vcv[4,3],
                            0,0,0,0,vcv[1,4],vcv[2,4],vcv[3,4],vcv[4,4]),
                          nrow=8,ncol=8)

        A <- matrix(c(1,0,0,0,0,0,
                      0,1,0,0,0,0,
                      0,0,1,0,0,0,
                      0,0,0,1,0,0,
                      1,1,0,0,0,0,
                      0,0,1,1,0,0,
                      0,0,0,0,1,0,
                      0,0,0,0,0,1),
                    nrow=6, ncol=8)

        sigma <- A%*%sigma88%*%t(A)

        return(sigma)
    }
}