compute_loglikeimp: Compute the log-likelihood of the observed data given PCA...

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

Description

The log-likelihood of the data for probabilistic PCA is known to be multivariate Gaussian. Using this, one can check the log-likelihood value of the observed data values given the parameter estimates from the PCA model. This can be useful to compare different models.

Usage

compute_loglikeimp(dat, A, S, covmat, meanvec, verbose = TRUE)

Arguments

dat	matrix – the data matrix with variables in rows and observations in columns.
A	matrix – estimated loadings matrix with observed variables in rows and latent variables in columns.
S	matrix – estimated factor scores matrix with latent variables in rows and observations in columns.
covmat	matrix – the estimated covariance matrix.
meanvec	numeric – the estimated mean vector.
verbose	logical – whether extra output should be displayed.

Value

the log-likelihood value

Examples

p <- 20
n <- 20
set.seed(10045)
  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(p*n), p, n)
  miss.inds <- sample(1:(p*n), round(p*n/10))
  X[miss.inds] <- NaN
  Xsaved <- X
  M <- !is.nan(X)
  X[X==0]      <- .Machine$double.eps
  X[is.nan(X)] <- 0
  
  notmiss <- which(X!=0, arr.ind = TRUE)
  IX      <- notmiss[,1]
  JX      <- notmiss[,2]
  
  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
  Mu <- rowSums(X) / Nobs_i
  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)
  # initialised model log-likelihood
  compute_loglikeimp(dat=Xsaved, A=A, S=S, covmat=tcrossprod(A)+diag(p),
  meanvec=Mu, verbose=TRUE)
  # estimated model log-likelihood
  compute_loglikeimp(dat=Xsaved, A=ppcaOutput$W, S=t(ppcaOutput$scores), covmat=ppcaOutput$C,
  meanvec=ppcaOutput$m, verbose=TRUE)

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