Nothing
#' Linearized Bregman solver for composite conditionally likelihood of Gaussian Graphical
#' model with lasso penalty.
#'
#' Solver for the entire solution path of coefficients.
#'
#' The data matrix X is assumed to follow the Gaussian Graohical model which is described as following:\cr
#' \deqn{X \sim N(\mu, \Theta^{-1})}\cr
#' where \eqn{\Theta} is sparse p-by-p symmetric matrix. Then conditional on \eqn{x_{-j}}\cr
#' \deqn{x_j \sim N(\mu_j - \sum_{k\neq j}\Theta_{jk}/\Theta_{jj}(x_k-\mu_k),1/\Theta_{jj}) }\cr
#' then the composite conditional likelihood is like this:\cr
#' \deqn{- \sum_{j} condloglik(X_j | X_{-j})}\cr
#' or in detail:\cr
#' \deqn{\sum_{j} \Theta_{j}^TS\Theta_{j}/2\Theta_{jj} - ln(\Theta_{jj})/2}\cr
#' where \eqn{S} is covariance matrix of data. It is easy to prove that this loss function
#' is convex.
#'
#'
#' @param X An n-by-p matrix of variables.
#' @param kappa The damping factor of the Linearized Bregman Algorithm that is
#' defined in the reference paper. See details.
#' @param alpha Parameter in Linearized Bregman algorithm which controls the
#' step-length of the discretized solver for the Bregman Inverse Scale Space.
#' See details.
#' @param S The covariance matrix can be provided directly if data matrix X is missing.
#' @param c Normalized step-length. If alpha is missing, alpha is automatically generated by
#' \code{alpha=c*n/(kappa*||X^T*X||_2)}. Default is 2. It should be in (0,4).
#' If beyond this range the path may be oscillated at large t values.
#' @param tlist Parameters t along the path.
#' @param nt Number of t. Used only if tlist is missing. Default is 100.
#' @param trate tmax/tmin. Used only if tlist is missing. Default is 100.
#' @param print If TRUE, the percentage of finished computation is printed.
#' @return A "ggm" class object is returned. The list contains the call,
#' the path, value for alpha, kappa, t.
#' @author Jiechao Xiong
#' @keywords regression
#' @examples
#'
#' library(MASS)
#' p = 20
#' Omega = diag(1,p,p)
#' Omega[0:(p-2)*(p+1)+2] = 1/3
#' Omega[1:(p-1)*(p+1)] = 1/3
#' S = solve(Omega)
#' X = mvrnorm(n=500,rep(0,p),S)
#' obj = ggm(X,10,trate=10)
#' obj$path[,,50]
ggm <- function(X, kappa, alpha, S=NA, c = 2, tlist,nt = 100,trate = 100,print=FALSE)
{
if (!missing(X)){
n <- nrow(X)
S <- cov(X)*(1-1/n)
}else if(any(is.na(S))){stop('A data matrix or covariance matrix must be provided!')}
p <- nrow(S)
if (missing(tlist)) tlist<-rep(-1.0,nt)
else nt<-length(tlist)
if (missing(alpha)){
sigma <- norm(S,"2")
alpha <- c/kappa/sigma^2
}
result_r <- vector(length = nt*p^2)
solution <- .C("ggm_C",
as.numeric(S),
as.integer(p),
as.numeric(kappa),
as.numeric(alpha),
as.numeric(result_r),
as.numeric(tlist),
as.integer(nt),
as.numeric(trate),
as.integer(print!=0)
)
path <- sapply(0:(nt- 1), function(x)
matrix(solution[[5]][(1+x*p^2):((x+1)*p^2)], p, p),simplify = "array")
object <- list(call = call,family="ggm", kappa = kappa, alpha = alpha, path = path, nt = nt,t=solution[[6]])
class(object) <- "ggm"
object
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.