pca_updates: PPCA updates

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

Description

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 pca_full.

Usage

pca_updates(
  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
)

Arguments

X	matrix – data matrix with variables in rows and observations in columns.
V	numeric – scalar value corresponding to the initialised variance of the error parameter.
A	matrix – initialised loadings matrix with observed variables in rows and latent variables in columns
Av	array – initialised covariance matrices of the rows of A.
Va	numeric – the hyperparameter of the prior variance of the rows of A.
S	matrix – initialised factor scores matrix with latent variables in rows and observations in columns.
Sv	array – initialised covariance matrices of the rows of S.
Mu	numeric – vector corresponding to the initialised mean of the observed variables.
Muv	numeric – the initialisation of the prior variance of Mu.
Vmu	numeric – the hyperparameter of the prior variance of Mu
hpVa	numeric – hyperparameter for the prior of the variance of Vmu and Va.
hpVb	numeric – hyperparameter for the prior of the variance of Vmu and Va.
hpV	numeric – hyperparameter for the prior of V.
ndata	numeric – number of observed values.
Nobs_i	numeric – number of observed values in each row of the data.
Isv	numeric – indices j for Svj that are identical. Not currently used.
M	matrix – logical values indicating which elements of the data are observed and missing.
IX	numeric – row indices of missing values
JX	numeric – column indices of missing values
rms	numeric – scalar indicating the initial rms
errMx	matrix – initial error matrix whose elements correspond to difference between the observed data and its model prediction.
bias	logical – value indicating whether the mean vector should be estimated or not.
rotate2pca	logical – value indicating whether to rotate the pca solution during learning.
niter_broadprior	numeric – number of iterations before the prior parameters begin to be updated.
use_prior	logical – whether or not a prior is assumed for the model parameters.
use_postvar	logical – whether the posterior variance should be computed and taken into account.
maxiters	numeric – the maximum number of iterations to be completed.
verbose	logical – whether extra output, such as the iteration number and cost function value, should be displayed.

Value

A list of 6 elements:

scores

matrix – the estimated scores.

m

numeric – the estimated mean vector.

ss

numeric – the estimated model variance.

W

matrix – the estimated loadings.

C

matrix – the estimated covariance matrix.

numIter

numeric – the number of iterations.

References

Ilin, A. and Raiko, T., 2010. link.

See Also

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)
  

HGray384/pcaNet documentation built on Nov. 14, 2020, 11:11 a.m.