patternCor3 <- function(A=diag(2), a=1:3){
# nk x nk covariance & correlation matrices
# k = length(a); abs(a) <= min(diag(A)
minV.A <- min(diag(A))
if(any(adj.a <- abs(a)>minV.A)){
warning("abs(a) too large; can't exceed",
" min(diag(A)) = ", minV.A,
"; forced into that range.")
a[adj.a] <- sign(a[adj.a])*minV.A
}
Aa <- kronecker(diag(3), A)
n <- dim(A)[1]
i1 <- n+1:n
i2 <- n+i1
diag(Aa[1:n, i1]) <- a[1]
diag(Aa[i1, 1:n]) <- a[1]
diag(Aa[1:n, i2]) <- a[2]
diag(Aa[i2, 1:n]) <- a[2]
diag(Aa[i1, i2]) <- a[3]
diag(Aa[i2, i1]) <- a[3]
s.A <- sqrt(diag(A))
r.Aa <- (Aa/outer(rep(s.A,3), rep(s.A,3)))
eig.Aa <- eigen(Aa)
list(Aa=Aa, corr.Aa=r.Aa, eigen.Aa=eig.Aa)
}If this works, all(eigen.Aa$values>=0). Thus, you can add a test for this and have something close to what you want. You could add an objective function that includes these eigenvalues with, say, minimum adjustment of "a" and feed it to "optim" and let "optim" find a solution that is "closest" in whatever sense you think is useful.
hope this helps. spencer graves
Gabor Grothendieck wrote:
Siew Leng TENG <siewlengteng <at> yahoo.com> writes:
: : Hi,
: : I would like to generate a correlation matrix with a
: particular structure. For example, a 3n x 3n matrix :
: A_(nxn) aI_(nxn) bI_(nxn)
: aI_(nxn) A_(nxn) cI_(nxn)
: aI_(nxn) cI_(nxn) A_(nxn)
: : where
: - A_(nxn) is a *specified* symmetric, positive
: definite nxn matrix.
: - I_(nxn) is an identity matrix of order n
: - a, b, c are (any) real numbers
: : Many attempts have been unsuccessful because a
: resulting matrix with any a, b, c may not be a
: positive definite one, and hence cannot qualify as a
: correlation matrix. Trying to first generate a
: covariance matrix however, does not guarantee a
: corresponding correlation matrix with the above
: structure.
: : My larger purpose is to use this correlation matrix to
: generate multivariate normal observations from the
: corresponding covariance matrix (derived via cholesky
: decomposition of the cor matrix).
This can be formulated a semidefinite programming problem.
I don't think R has any packages that do that but a google
search for "semidefinite programming" will find more info and some free non-R software which you could consider interfacing to R.
______________________________________________
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
-- Spencer Graves, PhD, Senior Development Engineer O: (408)938-4420; mobile: (408)655-4567
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
