R/RcppExports.R

In HGray384/pcaNet: Probabilistic principal components analysis - covariance estimation and network reconstruction

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' Bayesian PCA updates
#' 
#' Perform parameter updates for PPCA using the Variational Bayes framework
#' from Oba (2003). Not recommended to use standalone, rather it is called from within
#' \code{\link{bpcapM}} and its wrapper \code{\link{pcapM}}.
#' 
#' @param myMat \code{matrix} -- data matrix with observations in rows and 
#' variables in columns. (Note that this is the transpose of \code{X} in
#' \code{\link{pca_full}}.)
#' @param covy \code{matrix} -- the unbiased sample covariance of the data matrix.
#' @param N \code{numeric} -- the number of observations.
#' @param D \code{numeric} -- the number of variables.
#' @param hidden \code{numeric} -- indices of missing values in \code{1:length(myMat)}.
#' @param numberOfNonNAvaluesInEachCol \code{numeric} -- number of observed values in each 
#' column of the data (i.e. variables).
#' @param nomissIndex \code{numeric} -- indices of rows (observations) without any
#' missing values.
#' @param missIndex \code{numeric} -- indices of rows (observations) with missing
#' values.
#' @param nMissing \code{numeric} -- total number of missing values.
#' @param nPcs \code{numeric} -- number of components/latent variables to use.
#' @param threshold \code{numeric} -- threshold for convergence, applied to the precision
#' parameter \code{tau}. Updates for which the change in \code{tau} are below this threshold
#' value stop the algorithm.
#' @param maxIterations \code{numeric} -- the maximum number of iterations to be completed.
#' @return {A \code{list} of 5 elements:
#' \describe{
#' \item{W}{\code{matrix} -- the estimated loadings.}
#' \item{ss}{\code{numeric} -- the estimated model variance.}
#' \item{C}{\code{matrix} -- the estimated covariance matrix.}
#' \item{scores}{\code{matrix} -- the estimated scores.}
#' \item{m}{\code{numeric} -- the estimated mean vector.}
#' }}
#' @seealso \code{\link{bpcapM}}, \code{\link{pcapM}}
#' @examples
#' set.seed(102)
#'   N <- 20
#'   D <- 20
#'   nPcs <- 2
#'   maxIterations <- 1000
#'   X <- matrix(rnorm(50), D, N)
#'   X <- scale(X, center=TRUE, scale=FALSE) # mean 0
#'   covX <- cov(X)
#'   IX <- sample(1:D, 10)
#'   JX <- sample(1:N, 10)
#'   nMissing <- length(IX)+length(JX)
#'   X[JX, IX] <- 0
#'   hidden <- which(X==0)
#'   numberOfNonNAvaluesInEachCol <- colSums(X!=0)
#'   nomissIndex <- which(rowSums(X!=0)==N)
#'   missIndex <- which(rowSums(X!=0)!=N)
#'   threshold <- 1e-4
#'   bpcaNetOutput <- bpcaNet(myMat=X, covy=covX, N=N, D=D, hidden=hidden,
#'     numberOfNonNAvaluesInEachCol=numberOfNonNAvaluesInEachCol,
#'     nomissIndex=nomissIndex, missIndex=missIndex, nMissing=nMissing,
#'     nPcs=nPcs, threshold=threshold, maxIterations=maxIterations)
#' @references Oba, S., Sato, M.A., Takemasa, I., Monden, M.,
#'  Matsubara, K.I. and Ishii, S., 2003.
#'  \href{https://doi.org/10.1093/bioinformatics/btg287}{doi}.
#'  
#'  Stacklies, W., Redestig, H., Scholz, M., Walther, D. and 
#'  Selbig, J., 2007.
#'  \href{https://doi.org/10.1093/bioinformatics/btm069}{doi}.
bpcaNet <- function(myMat, covy, N, D, hidden, numberOfNonNAvaluesInEachCol, nomissIndex, missIndex, nMissing, nPcs = 2L, threshold = 1e-4, maxIterations = 200L) {
    .Call('_pcaNet_bpcaNet', PACKAGE = 'pcaNet', myMat, covy, N, D, hidden, numberOfNonNAvaluesInEachCol, nomissIndex, missIndex, nMissing, nPcs, threshold, maxIterations)
}

#' PPCA updates
#' 
#' Perform the parameter updates for PPCA using either Expectation-
#' Maximisation or Variational Bayes as in Ilin and Raiko (2010).
#' Recommended to not use standalone, rather this function is called
#' within \code{pca_full}.
#' 
#' @param X \code{matrix} -- data matrix with variables in rows and 
#' observations in columns.
#' @param V \code{numeric} -- scalar value corresponding to the
#' initialised variance of the error parameter.
#' @param A \code{matrix} -- initialised loadings matrix with observed
#' variables in rows and latent variables in columns
#' @param Av \code{array} -- initialised covariance matrices of the rows of
#' \code{A}.
#' @param Va \code{numeric} -- the hyperparameter of the prior variance of
#' the rows of \code{A}.
#' @param S \code{matrix} -- initialised factor scores matrix with latent
#' variables in rows and observations in columns.
#' @param Sv \code{array} -- initialised covariance matrices of the rows of 
#' \code{S}.
#' @param Mu \code{numeric} -- vector corresponding to the initialised mean 
#' of the observed variables.
#' @param Muv \code{numeric} -- the initialisation of the prior variance
#' of \code{Mu}.
#' @param Vmu \code{numeric} -- the hyperparameter of the prior variance 
#' of \code{Mu}
#' @param hpVa \code{numeric} -- hyperparameter for the prior of the variance
#' of \code{Vmu} and \code{Va}.
#' @param hpVb \code{numeric} -- hyperparameter for the prior of the variance
#' of \code{Vmu} and \code{Va}.
#' @param hpV \code{numeric} -- hyperparameter for the prior of
#' V.
#' @param ndata \code{numeric} -- number of observed values.
#' @param Nobs_i \code{numeric} -- number of observed values in each row
#' of the data.
#' @param Isv \code{numeric} -- indices j for Sv{j} that are identical. Not
#' currently used.
#' @param M \code{matrix} -- logical values indicating which elements of the
#' data are observed and missing.
#' @param IX \code{numeric} -- row indices of missing values
#' @param JX \code{numeric} -- column indices of missing values
#' @param rms \code{numeric} -- scalar indicating the initial rms
#' @param errMx \code{matrix} -- initial error matrix whose elements
#' correspond to difference between the observed data and its model 
#' prediction.
#' @param bias \code{logical} -- value indicating whether the mean
#' vector should be estimated or not.
#' @param rotate2pca \code{logical} -- value indicating whether to rotate
#' the pca solution during learning.
#' @param niter_broadprior \code{numeric} -- number of iterations before
#' the prior parameters begin to be updated. 
#' @param use_prior \code{logical} -- whether or not a prior is assumed
#' for the model parameters.
#' @param use_postvar \code{logical} -- whether the posterior variance
#' should be computed and taken into account.
#' @param maxiters \code{numeric} -- the maximum number of iterations to 
#' be completed.
#' @param verbose \code{logical} -- whether extra output, such as the
#' iteration number and cost function value, should be displayed.
#' 
#' @return {A \code{list} of 6 elements:
#' \describe{
#' \item{scores}{\code{matrix} -- the estimated scores.}
#' \item{m}{\code{numeric} -- the estimated mean vector.}
#' \item{ss}{\code{numeric} -- the estimated model variance.}
#' \item{W}{\code{matrix} -- the estimated loadings.}
#' \item{C}{\code{matrix} -- the estimated covariance matrix.}
#' \item{numIter}{\code{numeric} -- the number of iterations.}
#' }}
#' @seealso \code{\link{pca_full}}
#' @examples
#'   set.seed(102)
#'   n <- 20
#'   p <- 20
#'   verbose <- 1
#'   bias <- 1
#'   rotate2pca <- 1
#'   ncomp <- 2
#'   maxiters <- 1000
#'   opts <- list(init='random',
#'     maxiters=as.numeric(1000),
#'     niter_broadprior=as.numeric(100),
#'     earlystop=as.numeric(0)
#'   )
#'   use_prior = 1
#'   use_postvar = 1
#'   X <- matrix(rnorm(50), p, n)
#'   X <- t(scale(t(X), center=TRUE, scale=FALSE))
#'   IX <- sample(1:p, 10)
#'   JX <- sample(1:n, 10)
#'   X[IX, JX] <- 0
#'   M <- X!=0
#'   Nobs_i = rowSums(M)
#'   ndata   <- length(IX)
#'   # C++ indexing
#'   IX <- IX -1 
#'   JX <- JX -1
#'   
#'   initialisedParms <- initParms(p, n, ncomp, verbose = verbose)
#'   A   <- initialisedParms$A
#'   S   <- initialisedParms$S
#'   Mu  <- initialisedParms$Mu
#'   V   <- initialisedParms$V
#'   Av  <- initialisedParms$Av
#'   Sv  <- initialisedParms$Sv
#'   Muv <- initialisedParms$Muv
#'   Va  <- 1000*rep(1,ncomp)
#'   Vmu <- 1000
#'   if (is.null(Mu)){
#'     if (bias){
#'       Mu <- rowSums(X) / Nobs_i
#'   }else{
#'     Mu = rep(0, p)
#'     }
#'   }
#'   computedRMS <- compute_rms(X, A, S, M, ndata, verbose = verbose)
#'   errMx       <- computedRMS$errMx
#'   rms         <- computedRMS$rms
#'   hpVa <- 0.001
#'   hpVb <- 0.001
#'   hpV  <- 0.001
#'   Isv <- rep(0, 2)
#'   # data centering
#'   X <- subtractMu(Mu, X, M, p, n, bias, verbose = verbose) 
#'   ppcaOutput <- pca_updates(X=X, V=V, A=A, Va=Va, Av = Av, S = S, Sv = Sv, 
#'   Mu = Mu, Muv = Muv, Vmu = Vmu,
#'   hpVa = hpVa, hpVb = hpVb, hpV = hpV, ndata = ndata, Nobs_i = Nobs_i,
#'   Isv = Isv, M = M, IX = IX, JX = JX, rms = rms, errMx = errMx, 
#'   bias = bias, rotate2pca = rotate2pca, niter_broadprior = opts$niter_broadprior, 
#'   use_prior = use_prior, use_postvar = use_postvar,
#'   maxiters = maxiters, verbose = verbose)
#'   
#' @references Ilin, A. and Raiko, T., 2010.
#' \href{http://www.jmlr.org/papers/v11/ilin10a.html}{link}.
pca_updates <- function(X, V, A, Av, Va, S, Sv, Mu, Muv, Vmu, hpVa, hpVb, hpV, ndata, Nobs_i, Isv, M, IX, JX, rms, errMx, bias = 1L, rotate2pca = 1L, niter_broadprior = 100L, use_prior = 1L, use_postvar = 1L, maxiters = 1000L, verbose = 1L) {
    .Call('_pcaNet_pca_updates', PACKAGE = 'pcaNet', X, V, A, Av, Va, S, Sv, Mu, Muv, Vmu, hpVa, hpVb, hpV, ndata, Nobs_i, Isv, M, IX, JX, rms, errMx, bias, rotate2pca, niter_broadprior, use_prior, use_postvar, maxiters, verbose)
}

#' Probabilistic PCA updates
#' 
#' Perform parameter updates for PPCA using the Expectation-Maximisation framework
#' from Porta (2005) and also in the R-package \code{\link{pcaMethods}} (Stacklies, 2007).
#' Not recommended to use standalone, rather it is called from within
#' \code{\link{ppcapM}} and its wrapper \code{\link{pcapM}}.
#' 
#' @param myMat \code{matrix} -- data matrix with observations in rows and 
#' variables in columns. (Note that this is the transpose of \code{X} in
#' \code{\link{pca_full}}.)
#' @param N \code{numeric} -- the number of observations.
#' @param D \code{numeric} -- the number of variables.
#' @param W \code{matrix} -- initialised loadings matrix with observed
#' variables in rows and latent variables in columns.
#' @param hidden \code{numeric} -- indices of missing values in \code{1:length(myMat)}.
#' @param nMissing \code{numeric} -- total number of missing values.
#' @param nPcs \code{numeric} -- number of components/latent variables to use.
#' @param threshold \code{numeric} -- threshold for convergence, applied to the precision
#' parameter \code{tau}. Updates for which the change in \code{tau} are below this threshold
#' value stop the algorithm.
#' @param maxIterations \code{numeric} -- the maximum number of iterations to be completed.
#' 
#' @return {A \code{list} of 4 elements:
#' \describe{
#' \item{W}{\code{matrix} -- the estimated loadings.}
#' \item{ss}{\code{numeric} -- the estimated model variance.}
#' \item{C}{\code{matrix} -- the estimated covariance matrix.}
#' \item{myMat}{\code{matrix} -- the data matrix with missing values
#' replaced by their estimated projections.}
#' }}
#' @seealso \code{\link{ppcapM}}, \code{\link{pcapM}}
#' 
#' @references Porta, J.M., Verbeek, J.J. and Kroese, B.J., 2005.
#'  \href{https://hal.inria.fr/inria-00321476/en}{link}
#'  
#'  Stacklies, W., Redestig, H., Scholz, M., Walther, D. and 
#'  Selbig, J., 2007. \href{https://doi.org/10.1093/bioinformatics/btm069}{doi}.
#'  
#' @examples
#' set.seed(102)
#' N <- 20
#' D <- 20
#' nPcs <- 2
#' maxIterations <- 1000
#' X <- matrix(rnorm(50), D, N)
#' X <- scale(X, center=TRUE, scale=FALSE) # mean 0
#' covX <- cov(X)
#' IX <- sample(1:D, 10)
#' JX <- sample(1:N, 10)
#' nMissing <- length(IX)+length(JX)
#' X[JX, IX] <- 0
#' hidden <- which(X==0)
#' threshold <- 1e-4
#' r <- sample(N)
#' W <- t(X[r[1:nPcs], ,drop = FALSE])
#' W <- matrix(rnorm(W), nrow(W), ncol(W), dimnames = labels(W) )
#' ppcaNetOutput <- ppcaNet(X, N, D, W, hidden, nMissing, nPcs, threshold, maxIterations)
#' 
ppcaNet <- function(myMat, N, D, W, hidden, nMissing, nPcs = 2L, threshold = 1e-5, maxIterations = 1000L) {
    .Call('_pcaNet_ppcaNet', PACKAGE = 'pcaNet', myMat, N, D, W, hidden, nMissing, nPcs, threshold, maxIterations)
}