# 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

 1 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

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 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 Oct. 5, 2019, 3:36 p.m.