R/sigmafunc.R

SigmaGen <- function(vcv, ncauses){
    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)
}
kkholst/mcif documentation built on May 20, 2019, 10:47 a.m.