Nothing
#
# Implementation of Bayesian models for canonical correlation analysis (CCA),
# inter-battery factor analysis (BIBFA) and group factor analysis (GFA).
#
GFAexperiment <- function(Y,K,opts,Nrep=10) {
#
# A wrapper for running the GFA model Nrep times
# and choosing the final model based on the best
# lower bound. This is the recommended way of applying
# the algorithm.
#
# See GFA() for description of the inupts.
#
if(Nrep==1) {
return(GFA(Y,K,opts))
}
lb <- vector() # Lower bounds
models <- vector("list")
for(rep in 1:Nrep){
model <- GFA(Y,K,opts)
models[[rep]] <- model
lb[rep] <- tail(model$cost,1)
if(opts$verbose==1) {
print(paste("Run ",rep,"/",Nrep,": ",length(model$cost)," iterations",
" with final cost ",lb[rep],sep=""))
}
}
k <- which.max(lb)
model <- models[[k]]
return(model)
}
CCAexperiment <- function(Y,K,opts,Nrep=10) {
#
# A wrapper for running the CCA/BIBFA model Nrep times
# and choosing the final model based on the best
# lower bound. This is the recommended way of applying
# the algorithm.
#
# See CCA() for description of the inputs.
#
if(length(Y)!=2) {
print("Canonical correlation analysis (CCA) and Bayesian inter-battery
factor analysis (BIBFA) can only be computed for two co-occurring
data matrices. For multiple co-occurring data sources, use GFA()
and GFAexperiment() instead.")
} else {
return(GFAexperiment(Y,K,opts,Nrep))
}
}
CCA <- function(Y,K,opts) {
#
# A wrapper function for computing Bayesian CCA and BIBFA.
# Calls the same function that is used for GFA.
#
if(length(Y)!=2) {
print("Canonical correlation analysis (CCA) and Bayesian inter-battery
factor analysis (BIBFA) can only be computed for two co-occurring
data matrices. For multiple co-occurring data sources, use GFA()
and GFAexperiment() instead.")
} else {
return(GFA(Y,K,opts))
}
}
GFA <- function(Y,K,opts) {
#
# The main function for Bayesian group factor analysis
#
# Inputs:
# Y : List of M data matrices. Y[[m]] is a matrix with
# N rows (samples) and D_m columns (features). The
# samples need to be co-occurring.
# NOTE: All of these should be centered, so that the mean
# of each feature is zero
# NOTE: The algorithm is roughly invariant to the scale
# of the data, but extreme values should be avoided.
# Data with roughly unit variance or similar scale
# is recommended.
# K : The number of components
# opts : List of options (see function getDefaultOpts())
#
# Output:
# The trained model, which is a list that contains the following elements:
# Z : The mean of the latent variables; N times K matrix
# covZ : The covariance of the latent variables; K times K matrix
# ZZ : The second moments ZZ^T; K times K matrix
#
# W : List of the mean projections; D_i times K matrices
# covW : List of the covariances of the projections; D_i times D_i matrices
# WW : List of the second moments WW^T; K times K matrices
#
# tau : The mean precisions (inverse variance, so 1/tau gives the
# variances denoted by sigma in the paper); M-element vector
#
# alpha: The mean precisions of the projection weights, the
# variances of the ARD prior; M times K matrix
#
# U,V,u.mu,v.mu: The low-rank factorization of alpha.
#
# cost : Vector collecting the variational lower bounds for each
# iteration
# D : Data dimensionalities; M-element vector
# datavar : The total variance in the data sets, needed for
# GFAtrim()
# addednoise: The level of extra noise as in opts$addednoise
#
# Check that data is centered
if(!all(unlist(lapply(Y,colMeans)) < 1e-7) & opts$verbose == 2) {
print("Warning: The data does not have zero mean.")
}
# Check the number of views
if(length(Y)==1 & opts$verbose == 2){
print("Warning: The number of data sets must be larger than 1.")
}
#
# Store dimensionalities of data sets
#
M <- length(Y) # The number of views; M=2 corresponds to BIBFA
D <- unlist(lapply(Y,ncol)) # Collect the number of features in vector D
Ds <- sum(D) # Total number of features
N <- nrow(Y[[1]]) # The number of samples
datavar <- vector() # The total variance of the data, needed
for(m in 1:M) # for scaling in the initialization
datavar[m]=sum(apply(Y[[m]],2,var)) # and for GFAtrim()
R <- opts$R # Model rank R: either "full" or an integer from 0 to min(M,K)
if(R >= min(c(M,K)) && is.double(R)) { # if R >= min(M,K) convert to full
if(opts$verbose==2)
print("The rank corresponds to full rank solution.")
R <- "full"
}
if(R != "full") {
if(opts$verbose==2)
print("NOTE: Optimization of the rotation is not supported for low rank model.")
opts$rotate <- FALSE
}
# Some constants for speeding up the computation
const <- - N*Ds/2*log(2*pi) # Constant factors for the lower bound
Yconst <- unlist(lapply(Y,function(x){sum(x^2)}))
id <- rep(1,K) # Vector of ones for fast matrix calculations
alpha_0 <- opts$prior.alpha_0 # Easier access for hyperprior values
beta_0 <- opts$prior.beta_0
alpha_0t <- opts$prior.alpha_0t
beta_0t <- opts$prior.beta_0t
#
# Initialize the model randomly; other initializations could
# be done, but overdispersed random initialization is quite good.
#
# Latent variables
Z <- matrix(rnorm(N*K,0,1),N,K) # The mean
covZ <- diag(1,K) # The covariance
ZZ <- covZ + covZ*N # The second moments
# ARD and noise parameters
alpha <- matrix(1,M,K) # The mean of the ARD precisions
logalpha <- matrix(1,M,K) # The mean of <\log alpha >
if(R=="full") {
b_ard <- matrix(1,M,K) # The parameters of the Gamma distribution
a_ard <- alpha_0 + D/2 # for ARD precisions
digammaa_ard <- digamma(a_ard) # Constant for the lower bound
}
tau <- rep(opts$init.tau,M) # The mean noise precisions
a_tau <- alpha_0t + N*D/2 # The parameters of the Gamma distribution
b_tau <- rep(0,M) # for the noise precisions
digammaa_tau <- digamma(a_tau) # Constants needed for computing the lower bound
lgammaa_tau <- -sum(lgamma(a_tau))
lb.pt.const <- -M*lgamma(alpha_0t) + M*alpha_0t*log(beta_0t)
# Alpha needs to be initialized to match the data scale
for(m in 1:M) {
alpha[m,] <- K*D[m]/(datavar[m]-1/tau[m])
}
# The projections
# No need to initialize the projections randomly, since their updating
# step is the first one; just define the variables here
W <- vector("list",length=M) # The means
if(!opts$low.mem)
covW <- vector("list",length=M) # The covariances
else
covW <- diag(1,K)
WW <- vector("list",length=M) # The second moments
for(m in 1:M) {
W[[m]] <- matrix(0,D[m],K)
if(!opts$low.mem){
covW[[m]] <- diag(1,K)
WW[[m]] <- crossprod(W[[m]]) + covW[[m]]*D[m]
} else {
WW[[m]] <- crossprod(W[[m]]) + covW*D[m]
}
}
# Rotation parameters (full-rank only)
if(opts$rotate) {
Rot <- diag(K) # The rotation matrix
RotInv <- diag(K) # Its inverse
r <- as.vector(Rot) # Vectorized version of R
# parameter list for the optimization function (see ?optim)
par <- list(K=K,D=D,Ds=Ds,N=N,WW=WW,ZZ=ZZ,M=M)
}
# Use R-rank factorization of alpha.
if(R!="full") {
U <- matrix(abs(rnorm(M*R)),M,R)
lu <- length(U)
u.mu <- rep(0,M)
V <- matrix(abs((rnorm(K*R))),K,R)
lv <- length(V)
v.mu <- rep(0,K)
lambda <- opts$lambda
x <- c(as.vector(U),as.vector(V),u.mu,v.mu)
x <- rnorm(length(x))/100
par.uv <- list(getu=1:lu,getv=(lu+1):(lu+lv),getumean=(lu+lv+1):(lu+lv+M),
getvmean=(lu+lv+M+1):length(x),M=M,K=K,R=R,D=D,lambda=lambda)
par.uv$w2 <- matrix(0,M,K)
}
cost <- vector() # For storing the lower bounds
#
# The main loop
#
for(iter in 1:opts$iter.max) {
# Check if some components need to be removed
keep <- which(colMeans(Z^2) > 1e-7)
if(length(keep)!=K && opts$dropK) {
K <- length(keep)
if(K==0)
stop("Shut down all components, no structure found in the data.")
id <- rep(1,K)
Z <- Z[,keep,drop=F]
covZ <- covZ[keep,keep,drop=F]
ZZ <- ZZ[keep,keep,drop=F]
for(m in 1:M) {
W[[m]] <- W[[m]][,keep,drop=FALSE]
if(!opts$low.mem)
covW[[m]] <- covW[[m]][keep,keep,drop=F]
WW[[m]] <- WW[[m]][keep,keep,drop=F]
}
alpha <- alpha[,keep,drop=F]
logalpha <- logalpha[,keep,drop=F]
if(R!="full") {
V <- V[keep,,drop=F]
v.mu <- v.mu[keep]
x <- c(as.vector(U),as.vector(V),u.mu,v.mu)
lv <- length(V)
par.uv$K <- K
par.uv$getv<- (lu+1):(lu+lv)
par.uv$getumean <- (lu+lv+1):(lu+lv+M)
par.uv$getvmean <- (lu+lv+M+1):length(x)
par.uv$w2 <- matrix(0,M,K)
} else {
b_ard <- matrix(1,M,K) # The parameters of the Gamma distribution
}
if(opts$rotate) {
par$K <- K
}
}
#
# Update the projections
#
lb.qw <- rep(NA,M) # Computes also the determinant of covW needed for the
# lower bound
for(m in 1:M) {
# Efficient and robust way of computing
# solve(diag(alpha) + tau * ZZ^T)
tmp <- 1/sqrt(alpha[m,])
cho <- chol(outer(tmp,tmp)*ZZ + diag(1/tau[m],K))
det <- -2*sum(log(diag(cho))) - sum(log(alpha[m,])) - K*log(tau[m])
lb.qw[m] <- det
if(!opts$low.mem) {
covW[[m]] <- 1/tau[m] * outer(tmp,tmp) * chol2inv(cho)
W[[m]] <- crossprod(Y[[m]],Z)%*%covW[[m]]*tau[m]
WW[[m]] <- crossprod(W[[m]]) + covW[[m]]*D[m]
} else {
covW <- 1/tau[m] * outer(tmp,tmp) * chol2inv(cho)
W[[m]] <- crossprod(Y[[m]],Z)%*%covW*tau[m]
WW[[m]] <- crossprod(W[[m]]) + covW*D[m]
}
}
#
# Update the latent variables
#
# Efficient and robust way of computing
# solve(diag(1,K) + tau * WW^t)
covZ <- diag(1,K)
for(m in 1:M) {
covZ <- covZ + tau[m]*WW[[m]]
}
cho <- chol(covZ)
covZ <- chol2inv(cho)
det <- -2*sum(log(diag(cho)))
lb.qx <- det
Z <- Z*0
for(m in 1:M)
Z <- Z + Y[[m]]%*%W[[m]]*tau[m]
Z <- Z%*%covZ
ZZ <- crossprod(Z) + N*covZ
#
# Optimization of the rotation (only start after the first
# iteration)
#
if(R=="full" & opts$rotate & iter > 1) {
# Update the parameter list for the optimizer
par$WW <- WW
par$ZZ <- ZZ
# Always start from the identity matrix, i.e. no rotation
r <- as.vector(diag(K))
if(opts$opt.method == "BFGS") {
r.opt <- try(optim(r,E,gradE,par,method="BFGS",
control=list(reltol=opts$opt.bfgs.crit,
maxit=opts$opt.iter)), silent=TRUE)
}
if(opts$opt.method== "L-BFGS") {
r.opt <- try(optim(r,E,gradE,par,method="L-BFGS-B",
control=list(maxit=opts$opt.iter,
factr=opts$lbfgs.factr)), silent=TRUE)
}
# For large K the optimizer occasionally fails due to problems
# in the svd() routine
if(inherits(r.opt,"try-error")) {
print("Failure in optimizing the rotation. Turning the rotation off.")
opts$rotate <- FALSE
}else{
# Update the parameters involved in the rotation:
Rot <- matrix(r.opt$par,K)
eS <- svd(Rot)
det <- sum(log(eS$d))
RotInv <- tcrossprod( eS$v*outer(id, 1/eS$d), eS$u)
Z <- tcrossprod(Z,RotInv)
covZ <- tcrossprod( RotInv%*%covZ, RotInv)
ZZ <- crossprod(Z) + N*covZ
lb.qx <- lb.qx - 2*det
for(m in 1:M) {
if(!opts$low.mem) {
W[[m]] <- W[[m]]%*%Rot
covW[[m]] <- crossprod(Rot,covW[[m]])%*%Rot
WW[[m]] <- crossprod(W[[m]]) + covW[[m]]*D[m]
} else {
#covW[[m]] is not stored, so it needs to be computed before rotation
covW <- (WW[[m]] - crossprod(W[[m]]))/D[m]
W[[m]] <- W[[m]]%*%Rot
covW <- crossprod(Rot,covW)%*%Rot
WW[[m]] <- crossprod(W[[m]]) + covW*D[m]
}
lb.qw[m] <- lb.qw[m] + 2*det
}
}
}
#
# Update alpha, the ARD parameters
#
if(R == "full"){
for(m in 1:M){
tmp <- beta_0 + diag(WW[[m]])/2
alpha[m,] <- a_ard[m]/tmp
b_ard[m,] <- tmp
}
}else{
for(m in 1:M)
par.uv$w2[m,] <- diag(WW[[m]])
res <- try(optim(x,fn=Euv,gr=gradEuv,par.uv,method="L-BFGS-B",
lower=c(rep(-sqrt(500/R),M*R+K*R),rep(-50,M+K)),
upper=c(rep(sqrt(500/R),M*R+K*R), rep(50,M+K)),
control=list(fnscale=-1,maxit=opts$opt.iter,
factr=opts$lbfgs.factr)),silent=TRUE)
if(inherits(res,"try-error")){
cost[iter] <- NA
stop("Problems in optimization. Try a new initialization.")
# terminate the algorithm (next model to learn)
}
x <- res$par
U <- matrix(res$par[par.uv$getu],par.uv$M,par.uv$R)
V <- matrix(res$par[par.uv$getv],par.uv$K,par.uv$R)
u.mu <- res$par[par.uv$getumean]
v.mu <- res$par[par.uv$getvmean]
alpha <- exp(tcrossprod(U,V) + outer(u.mu,rep(1,K)) +outer(rep(1,M),v.mu))
}
#
# Update tau, the noise precisions
#
for(m in 1:M) {
b_tau[m] <- beta_0t + ( Yconst[m] + sum(WW[[m]]*ZZ)
- 2*sum(Z*(Y[[m]]%*%W[[m]])))/2
}
tau <- a_tau/b_tau
#
# Calculate the lower bound.
# Consists of calculating the likelihood term and
# KL-divergences between the factorization and the priors
#
# The likelihood
logtau <- digammaa_tau - log(b_tau)
if(R == "full") {
for(m in 1:M) {
logalpha[m,] <- digammaa_ard[m] - log(b_ard[m,])
}
} else {
logalpha <- log(alpha)
}
lb.p <- const + N*crossprod(D,logtau)/2 - crossprod(b_tau - beta_0t,tau)
lb <- lb.p
# E[ ln p(Z)] - E[ ln q(Z) ]
lb.px <- - sum(diag(ZZ))/2
lb.qx <- - N*lb.qx/2 - N*K/2
lb <- lb + lb.px - lb.qx
# E[ ln p(W)] - E[ ln q(W)]
if(R == "full") {
lb.pw <- 0
for(m in 1:M) {
lb.pw <- lb.pw +
D[m]/2*sum(logalpha[m,]) - sum(diag(WW[[m]])*alpha[m,])/2
}
} else {
lb.pw <- res$value
}
for(m in 1:M)
lb.qw[m] <- - D[m]*lb.qw[m]/2 -D[m]*K/2
lb <- lb + lb.pw - sum(lb.qw)
# E[ln p(alpha)] - E[ln q(alpha)]
if(R == "full") {
lb.pa <- M*K*( -lgamma(alpha_0) + alpha_0*log(beta_0) ) +
(alpha_0-1)*sum(logalpha) - beta_0*sum(alpha)
lb.qa <- -K*sum(lgamma(a_ard)) + sum(a_ard*rowSums( log(b_ard) )) +
sum((a_ard-1)*rowSums(logalpha)) - sum(b_ard*alpha)
lb <- lb + lb.pa - lb.qa
}
# E[ln p(tau)] - E[ln q(tau)]
lb.pt <- lb.pt.const + sum((alpha_0t-1)*logtau) - sum(beta_0t*tau)
lb.qt <- lgammaa_tau + crossprod(a_tau,log(b_tau)) +
crossprod((a_tau-1),logtau) - crossprod(b_tau,tau)
lb <- lb + lb.pt - lb.qt
# Store the cost function
cost[iter] <- lb
if(opts$verbose==2) {
print(paste("Iteration:",iter,"/ cost:",lb,"/ K",K))
}
# Convergence if the relative change in cost is small enough
if(iter>1){
diff <- cost[iter]-cost[iter-1]
if(abs(diff)/abs(cost[iter]) < opts$iter.crit | iter == opts$iter.max) {
break
}
}
} # the main loop of the algorithm ends
# Add a tiny amount of noise on top of the latent variables,
# to supress possible artificial structure in components that
# have effectively been turned off
Z <- Z + opts$addednoise*matrix(rnorm(N*K,0,1),N,K) %*% chol(covZ)
# return the output of the model as a list
if(R=="full") {
list(W=W,covW=covW,ZZ=ZZ,WW=WW,Z=Z,covZ=covZ,tau=tau,alpha=alpha,cost=cost,
D=D,K=K,addednoise=opts$addednoise,datavar=datavar,iter=iter,R=R)
} else {
list(W=W,covW=covW,ZZ=ZZ,WW=WW,Z=Z,covZ=covZ,tau=tau,alpha=alpha,cost=cost,
D=D,K=K,addednoise=opts$addednoise,datavar=datavar,iter=iter,R=R,
U=U,V=V,u.mu=u.mu,v.mu=v.mu)
}
}
E <- function(r,par) {
#
# Evaluates the (negative) cost function value wrt the transformation
# matrix R used in the generic optimization routine
#
R <- matrix(r,par$K)
eS <- svd(R)
val <- -sum(par$ZZ*tcrossprod(eS$u*outer(rep(1,par$K),1/eS$d)))/2
val <- val + (par$Ds-par$N)*sum(log(eS$d))
for(m in 1:par$M) {
val <- val - par$D[m]*sum( log( colSums(R*(par$WW[[m]]%*%R)) ))/2
}
return(-val)
}
gradE <- function(r,par) {
#
# Evaluates the (negative) gradient of the cost function E()
#
R <- matrix(r,par$K)
eS <- svd(R)
Rinv <- tcrossprod(eS$v*outer(rep(1,par$K),1/eS$d), eS$u)
gr <- as.vector(tcrossprod(tcrossprod(eS$u*outer(rep(1,par$K),1/eS$d^2),eS$u)%*%
par$ZZ + diag(par$Ds - par$N,par$K), Rinv) )
tmp1 <- par$WW[[1]]%*%R
tmp2 <- 1/colSums(R*tmp1)
tmp1 <- par$D[1]*as.vector( tmp1*outer(rep(1,par$K),tmp2) )
gr <- gr - tmp1
for(m in 2:par$M){
tmp1 <- par$WW[[m]]%*%R
tmp2 <- 1/colSums(R*tmp1)
tmp1 <- par$D[m]*as.vector( tmp1*outer(rep(1,par$K),tmp2) )
gr <- gr - tmp1
}
return(-gr)
}
Euv <- function(x,par){
#
# Evaluates the cost function value wrt the low-rank
# factorization of alpha used in the generic optimization routine
#
U <- matrix(x[par$getu],par$M,par$R)
V <- matrix(x[par$getv],par$K,par$R)
u.mu <- x[par$getumean]
v.mu <- x[par$getvmean]
logalpha <- tcrossprod(U,V) + outer(u.mu,rep(1,par$K)) + outer(rep(1,par$M),v.mu)
E <- sum(crossprod(par$D,logalpha))-sum(par$w2*exp(logalpha))
if(par$lambda!=0)
E <- E - par$lambda*(sum(V^2) + sum(U^2))
return(E/2)
}
gradEuv <- function(x,par){
#
# Evaluates the gradient of the cost function Euv()
#
U <- matrix(x[par$getu],par$M,par$R)
V <- matrix(x[par$getv],par$K,par$R)
u.mu <- x[par$getumean]
v.mu <- x[par$getvmean]
alphaiAlphaw2 <- outer(par$D,rep(1,par$K)) -
exp(tcrossprod(U,V) + outer(u.mu,rep(1,par$K)) + outer(rep(1,par$M),v.mu))*par$w2
gradU <- alphaiAlphaw2%*%V
gradV <- crossprod(alphaiAlphaw2,U)
if(par$lambda!=0) {
gradU <- gradU - par$lambda*2*U
gradV <- gradV - par$lambda*2*V
}
grad.umean <- rowSums(alphaiAlphaw2)
grad.vmean <- colSums(alphaiAlphaw2)
grad <- c(as.vector(gradU),as.vector(gradV),grad.umean,grad.vmean)
return(grad/2)
}
getDefaultOpts <- function(){
#
# A function for generating a default set of parameters.
#
# To run the algorithm with other values:
# opts <- getDefaultOpts()
# opts$opt.method <- "BFGS"
# model <- GFA(Y,K,opts)
#
# Rank of GFA's hierarhical low-rank ARD prior. Possible values are "full"
# and all integers, including zero.
#
R <- "full"
#
# Prior precision of zero-mean Gaussian prior for U and V of
# the low-rank ARD prior.
#
lambda <- 0.1
#
# Whether to use the rotation explained in the ICML'11 paper.
# Using the rotation is recommended, but for some data sets
# equally good solutions can be obtained without rotation and
# possibly faster.
# - TRUE|FALSE
#
rotate <- TRUE
#
# Parameters for controlling how the rotation is solved
# - opt.method chooses the optimization method and
# takes values "BFGS" or "L-BFGS". The former
# is typically faster but takes more memory, so the latter
# is the default choice. For small K may use BFGS instead.
# - opt.iter is the maximum number of iterations
# - lbfgs.factr is convergence criterion for L-BFGS; smaller
# values increase the accuracy (10^7 or 10^10 could be tried
# to speed things up)
# - bfgs.crit is convergence criterion for BFGS; smaller
# values increase the accuracy (10^-7 or 10^-3 could also be used)
#
opt.method <- "L-BFGS"
opt.iter <- 10^5
lbfgs.factr <- 10^10
bfgs.crit <- 10^-5
#
# Initial value for the noise precisions. Should be large enough
# so that the real structure is modeled with components
# instead of the noise parameters (see Luttinen&Ilin, 2010)
# Values: Positive numbers, but generally should use values well
# above 1
#
init.tau <- 10^3
#
# Parameters for controlling when the algorithm stops.
# It stops when the relative difference in the lower bound
# falls below iter.crit or iter.max iterations have been performed.
#
iter.crit <- 10^-6
iter.max <- 10^5
#
# Additive noise level for latent variables. The latent variables
# of inactive components (those with very large alpha) occasionally
# show some structure in the mean values, even though the distribution
# matches very accurately the prior N(0,I). This structure disappears
# is a tiny amount of random noise is added on top of the
# mean estimates. Setting the value to 0 will make the predictions
# deterministic
#
addednoise <- 1e-5
#
# Hyperparameters
# - alpha_0, beta_0 for the ARD precisions
# - alpha_0t, beta_0t for the residual noise predicions
#
prior.alpha_0 <- prior.beta_0 <- 1e-14
prior.alpha_0t <- prior.beta_0t <- 1e-14
# Two performace enhancing settings
# dropK=TRUE : matrix dimensions are reduced when components are shut off
# low.mem=TRUE: list covW is not stored (may cause minor numerical differences)
dropK <- TRUE
low.mem <- FALSE
#
# Verbosity level
# 0: Nothing
# 1: Final cost function value for each run of GFAexperiment()
# 2: Cost function values for each iteration
#
verbose <- 2
return(list(R=R,rotate=rotate, init.tau=init.tau, iter.crit=iter.crit,
iter.max=iter.max, opt.method=opt.method,
lbfgs.factr=lbfgs.factr, bfgs.crit=bfgs.crit, opt.iter=opt.iter,
addednoise=addednoise,
prior.alpha_0=prior.alpha_0,prior.beta_0=prior.beta_0,
prior.alpha_0t=prior.alpha_0t,prior.beta_0t=prior.beta_0t,
verbose=verbose,lambda=lambda,dropK=dropK,low.mem=low.mem))
}
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.