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)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.