Message-ID: <CANz9Z_+e1_wsPi3r833KY8QZAscyN3WC7V78hkhWdNHc+qvNgA@mail.gmail.com>
Date: 2011-10-19T03:31:32Z
From: Joshua Wiley
Subject: Sparse covariance estimation (via glasso) shrinking to a "nonzero" constant
In-Reply-To: <20111018192431.11742fguxhn9mask@mail.ucla.edu>
Hi Vivian,
This may be naive given the method (I am unfamiliar with glasso), but
what about simple subtraction? If it restricts to 0, you believe you
have .3, then just: obs - .3 and restrict to 0 again? Here is a
little example (assuming .3 correlation, but using glasso with the
covariance matrix so subtraction is marginally more complex):
dat <- lapply(list((dat <- prcomp(matrix(rnorm(1000), 200))$x), dat %*%
chol(matrix(c(1, rep(c(rep(.3, 5), 1), 4)), 5, 5))), cov)
i <- sqrt(1/diag(dat[[2]]))
v <- dat[[2]] - matrix(c(0, rep(c(rep(.3, 5), 0), 4)), 5, 5)/(i *
rep(i, each = 5))
require(glasso)
Reduce(`-`, res <- list(true = glasso(dat[[1]], .01)$w, fake =
glasso(v, .01)$w))
Again this may not make any sense at all given how the methods actually work.
Cheers,
Josh
On Tue, Oct 18, 2011 at 7:24 PM, Vivian Shih <vivs at ucla.edu> wrote:
> I've only been using R on and off for 9 months and started using the glasso
> package for sparse covariance estimation. I know the concept is to shrink
> some of the elements of the covariance matrix to zero.
>
> However, say I have a dataset that I know has some underlying "baseline"
> covariance/correlation (say, a value of 0.3), how can I change or
> incorporate that into to the code? Basically instead of "let's restrict some
> of these elements are 0 and others have different nonzero values", I'd like
> it to be "let's restrict some of these elements are 0.3 and others are
> different or higher/lower than this 0.3". Does that make any sense?
>
>
>
> Here is the glasso code if it helps:
>>
>> library(glasso)
>> glasso
>
> function (s, rho, zero = NULL, thr = 1e-04, maxit = 10000, approx = FALSE,
> ? ?penalize.diagonal = TRUE, start = c("cold", "warm"), w.init = NULL,
> ? ?wi.init = NULL, trace = FALSE)
> {
> ? ?n = nrow(s)
> ? ?BIG = 1e+10
> ? ?if (!is.matrix(rho) & length(rho) != 1 & length(rho) != nrow(s)) {
> ? ? ? ?stop("Wrong number of elements in rho")
> ? ?}
> ? ?if (length(rho) == 1 & rho == 0) {
> ? ? ? ?cat("With rho=0, there may be convergence problems if the input
> matrix s\n is not of full rank",
> ? ? ? ? ? ?fill = TRUE)
> ? ?}
> ? ?if (is.vector(rho)) {
> ? ? ? ?rho = matrix(sqrt(rho)) %*% sqrt(rho)
> ? ?}
> ? ?if (length(rho) == 1) {
> ? ? ? ?rho = matrix(rho, ncol = n, nrow = n)
> ? ?}
> ? ?if (!is.null(zero)) {
> ? ? ? ?if (!is.matrix(zero)) {
> ? ? ? ? ? ?zero = matrix(zero, nrow = TRUE)
> ? ? ? ?}
> ? ? ? ?for (k in 1:nrow(zero)) {
> ? ? ? ? ? ?i = zero[k, 1]
> ? ? ? ? ? ?j = zero[k, 2]
> ? ? ? ? ? ?rho[i, j] = BIG
> ? ? ? ? ? ?rho[j, i] = BIG
> ? ? ? ?}
> ? ?}
> ? ?start.type = match.arg(start)
> ? ?if (start.type == "cold") {
> ? ? ? ?is = 0
> ? ? ? ?ww = xx = matrix(0, nrow = n, ncol = n)
> ? ?}
> ? ?if (start.type == "warm") {
> ? ? ? ?is = 1
> ? ? ? ?if (is.null(w.init) | is.null(wi.init)) {
> ? ? ? ? ? ?stop("Warm start specified: w.init and wi.init must be non-null")
> ? ? ? ?}
> ? ? ? ?ww = w.init
> ? ? ? ?xx = wi.init
> ? ?}
> ? ?itrace = 1 * trace
> ? ?ipen = 1 * (penalize.diagonal)
> ? ?ia = 1 * approx
> ? ?mode(rho) = "double"
> ? ?mode(s) = "double"
> ? ?mode(ww) = "double"
> ? ?mode(xx) = "double"
> ? ?mode(n) = "integer"
> ? ?mode(maxit) = "integer"
> ? ?mode(ia) = "integer"
> ? ?mode(itrace) = "integer"
> ? ?mode(ipen) = "integer"
> ? ?mode(is) = "integer"
> ? ?mode(thr) = "double"
> ? ?junk <- .Fortran("glasso", n, s, rho, ia, is, itrace, ipen,
> ? ? ? ?thr, maxit = maxit, ww = ww, xx = xx, niter = integer(1),
> ? ? ? ?del = double(1), ierr = integer(1), PACKAGE = "glasso")
> ? ?ww = matrix(junk$ww, ncol = n)
> ? ?xx = matrix(junk$xx, ncol = n)
> ? ?if (junk$ierr != 0) {
> ? ? ? ?stop("memory allocation error")
> ? ?}
> ? ?critfun = function(Sigmahati, s, rho, penalize.diagonal = TRUE) {
> ? ? ? ?d = det(Sigmahati)
> ? ? ? ?temp = Sigmahati
> ? ? ? ?if (!penalize.diagonal) {
> ? ? ? ? ? ?diag(temp) = 0
> ? ? ? ?}
> ? ? ? ?val = -log(d) + sum(diag(s %*% Sigmahati)) + sum(abs(rho *
> ? ? ? ? ? ?temp))
> ? ? ? ?return(val)
> ? ?}
> ? ?crit = critfun(xx, s, rho, penalize.diagonal)
> ? ?return(list(w = ww, wi = xx, loglik = -(n/2) * crit, errflag = junk$ierr,
> ? ? ? ?approx = approx, del = junk$del, niter = junk$niter))
> }
>
>
>
> If anyone can give any suggestions, I'd greatly appreciate it. Thanks!
>
>
>
>
> Vivian
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
--
Joshua Wiley
Ph.D. Student, Health Psychology
Programmer Analyst II, ATS Statistical Consulting Group
University of California, Los Angeles
https://joshuawiley.com/