gapsRun: `gapsRun` calls the C++ MCMC code and performs Bayesian matrix factorization returning the two matrices that reconstruct the data matrix

Description

gapsRun calls the C++ MCMC code and performs Bayesian matrix factorization returning the two matrices that reconstruct the data matrix

Usage

gapsRun(D, S, ABins = data.frame(), PBins = data.frame(), nFactor = 7,
  simulation_id = "simulation", nEquil = 1000, nSample = 1000,
  nOutR = 1000, output_atomic = "FALSE", fixedBinProbs = "FALSE",
  fixedDomain = "N", sampleSnapshots = "TRUE", numSnapshots = 100,
  alphaA = 0.01, nMaxA = 1e+05, max_gibbmass_paraA = 100, alphaP = 0.01,
  nMaxP = 1e+05, max_gibbmass_paraP = 100)

Arguments

data matrix

uncertainty matrix (std devs for chi-squared of Log Likelihood)

ABins

a matrix of same size as A which gives relative probability of that element being non-zero (optional)

PBins

a matrix of same size as P which gives relative probability of that element being non-zero (optional)

nFactor

number of patterns (basis vectors, metagenes), which must be greater than or equal to the number of rows of FP (optional, defaults to 7)

simulation_id

name to attach to atoms files if created

nEquil

number of iterations for burn-in

nSample

number of iterations for sampling

nOutR

how often to print status into R by iterations

output_atomic

whether to write atom files (large)

fixedBinProbs

Boolean for using relative probabilities given in Abins and Pbins

fixedDomain

character to indicate whether A or P is domain for relative probabilities

sampleSnapshots

Boolean to indicate whether to capture individual samples from Markov chain during sampling

numSnapshots

the number of individual samples to capture

alphaA

sparsity parameter for A domain

nMaxA

PRESENTLY UNUSED, future = limit number of atoms

max_gibbmass_paraA

limit truncated normal to max size

alphaP

sparsity parameter for P domain

nMaxP

PRESENTLY UNUSED, future = limit number of atoms

max_gibbmass_paraP

limit truncated normal to max size

Value

A list containing:
AmeanSampled mean value of the amplitude matrix ${\bf{A}}$.
AsdSampled standard deviation of the amplitude matrix ${\bf{A}}$.
PmeanSampled mean value of the amplitude matrix ${\bf{P}}$.
PsdSampled standard deviation of the amplitude matrix ${\bf{P}}$.
atomsAEquilNumber of atoms in ${\bf{A}}$ during each iteration of the equilibration phase.
atomsASampNumber of atoms in ${\bf{A}}$ during each iteration of the sampling phase.
atomsPEquilNumber of atoms in ${\bf{P}}$ during each iteration of the equilibration phase.
atomsPSampNumber of atoms in ${\bf{P}}$ during each iteration of the sampling phase.
chiSqValuesValue of $chi^2$ at each step during equilibration and sampling.
meanChi2Value of $chi^2$ for Amean and Pmean.
ASnapshotsSamples of A matrices taken during sampling.
PSnapshotsSamples of P matrices taken during sampling.

Details

The decomposition in GAPS is achieved by finding amplitude and pattern matrices (${\bf{A}}$ and ${\bf{P}}$, respectively) for which $${\bf{D}} = {\bf{A}}{\bf{P}} + \Sigma$$, where $\Sigma$ is the matrix of uncertainties given by S. The matrices $\bf{A}$ and $\bf{P}$ are assumed to have the atomic prior described in Sibisi and Skilling (1997) and are found with MCMC sampling.

Examples

Run this code

## Load data
data('SimpSim')

## Run GAPS matrix decomposition
nIter <- 5000
results <- gapsRun(SimpSim.D, SimpSim.S, nFactor=3,
                   nEquil=nIter, nSample=nIter)

## Plot the results
plotGAPS(results$Amean, results$Pmean, 'GSFigs')