Compute the posterior joint probability distribution of the variates conditional on the given data, by means of Markov-chain Monte Carlo, using the package Nimble.
learn(
data,
metadata,
auxdata = NULL,
outputdir = NULL,
nsamples = 3600,
nchains = 8,
nsamplesperchain = 450,
parallel = TRUE,
seed = NULL,
cleanup = TRUE,
appendinfo = TRUE,
valueislearnt = TRUE,
subsampledata = NULL,
prior = missing(data) || is.null(data),
startupMCiterations = 3600,
minMCiterations = 0,
maxMCiterations = +Inf,
maxhours = +Inf,
ncheckpoints = 12,
maxrelMCSE = +Inf,
minESS = 450,
initES = 2,
thinning = NULL,
verbose = TRUE,
plottraces = !cleanup,
showKtraces = FALSE,
showAlphatraces = FALSE,
hyperparams = list(ncomponents = 64, minalpha = -4, maxalpha = 4, byalpha = 1, Rshapelo
= 0.5, Rshapehi = 0.5, Rvarm1 = 3^2, Cshapelo = 0.5, Cshapehi = 0.5, Cvarm1 = 3^2,
Dshapelo = 0.5, Dshapehi = 0.5, Dvarm1 = 3^2, Bshapelo = 1, Bshapehi = 1, Dthreshold
= 1, tscalefactor = 4.266, Oprior = "Hadamard", Nprior = "Hadamard", avoidzeroW =
NULL, initmethod = "datacentre", Qerror = pnorm(c(-1, 1)))
)A "learnt" object, or name of directory containing such an object and other output files, or NULL, depending on argument valueislearnt.
learn() saves several files in a directory. By default this output directory is a temporary directory within the one used by base::tempdir(), but an alternative one can be chosen with the argument outputdir =. The output directory contain several diagnostic files for the Monte Carlo computation; in particular:
MCtraces.pdf: shows several trace plots of the Monte Carlo sampling; the correspondin data are in the file MCtraces.rds.
plotsamples_learnt.pdf, plotquantiles_learnt.pdf: show the marginal posterior distributions of each individual variate, together with their variability (as samples or quantiles).
log-1.log, log-2.log, ... one for each parallel core; report the progress of each parallel Monte Carlo computation and notes about it.
rng_seed.rds: the state of the pseudorandom seed (see base::Random) when learn() was called.
metadata.csv: a copy of the metadata.
It is recommended that you give an explicit argument outputdir = and save the directory with the files above for future reference. In particular, the MCtraces.pdf plot and MCtraces.rds data can be useful to report Monte Carlo convergence in any work of yours that used Prova.
A dataset, given as a base::data.frame() or as a file path to a CSV file. If missing or NULL, then the prior probability distribution is calculated.
metadata about the dataset's variates, given either as a data frame or as a file path to a CSV file.
A larger dataset, given as a data frame or as a file path to a CSV file. Such a dataset would be too large to use in the Monte Carlo sampling, but can still be used to help estimate some hyperparameters.
NULL (default) or NA or character: path to folder where output information and diagnostics should be saved. If NULL, a directory is created in the temporary-directory space given by base::tempdir(). If NA, a directory is created in the current working directory given by base::getwd(). If character, this is taken to be the output directory; it should of course be writable by the user.
Integer, default 3600: number of desired, approximately independent Monte Carlo samples. If this argument is changed, the user is also required to explicitly give either nchains or nsamplesperchain, but not both; the remaining third argument is determined from \(\mathrm{nsamples} = \mathrm{nchains} \times \mathrm{nsamplesperchain}\).
Integer, default 8: number of Monte Carlo chains. If this argument is changed, the user is also required to explicitly give either nsamples or nsamplesperchain, but not both; the remaining third argument is determined from \(\mathrm{nsamples} = \mathrm{nchains} \times \mathrm{nsamplesperchain}\).
Integer, default 450: number of approximately independent Monte Carlo samples per chain. If this argument is changed, the user is also required to explicitly give either nsamples or nchains, but not both; the remaining third argument is determined from \(\mathrm{nsamples} = \mathrm{nchains} \times \mathrm{nsamplesperchain}\).
Logical or positive integer or cluster object. TRUE (default): use roughly half of available cores; FALSE (default): use serial computation; integer: use this many cores. It can also be a cluster object previously created with parallel::makeCluster(); in this case the parallel computation will use this object.
Integer or NULL (default): use this seed for the random number generator. If NULL, do not set the seed.
Logical, default TRUE: remove diagnostic files at the end of the computation?
Logical, default TRUE: append information about number of variates ('V'), number of data points ('D'), number of Monte Carlo samples ('S'), and timestamp, to the name of the output directory outputdir? The appended string has the format 'Vn_Dn_Sn_YYMMDDTHHMMSS'.
Logical or NULL: should the VALUE returned be the learnt object containing the results from the Monte Carlo computation? Default TRUE. If FALSE, then VALUE is the output directory name. If NULL, then VALUE is NULL.
Integer or NULL (default): if integer, use only that many datapoints from the original dataset in the data argument.
Logical: Calculate the prior distribution? Default is FALSE unless data argument is missing or NULL.
Integer, default 3600: number of initial Monte Carlo iterations.
Integer, default 0: minimum number of Monte Carlo iterations to be doneby a chain.
Integer, default Inf: Do at most this many Monte Carlo iterations per chain.
Numeric, default Inf: approximate time limit, in hours, for the Monte Carlo computation to last.
Integer or NULL, default 12: number of datapoints (per chain) to use for checking when the Monte Carlo computation should end. If NULL, this is equal to number of variates + 2. If Inf, use all datapoints.
Numeric positive, default +Inf: desired maximal relative Monte Carlo Standard Error of calculated probabilities with respect to their variability with new data. The default +Inf means that minESS is used instead. maxrelMCSE is related to minESS by \(\mathrm{maxrelMCSE} = 1/\sqrt{\mathrm{minESS} + \mathrm{initES}}\).
Numeric positive or NULL, default 450: desired minimal Monte Carlo Expected Sample Size. If NULL, it is equal to the final nsamplesperchain. minESS is related to maxrelMCSE by \(\mathrm{minESS} = 1/\mathrm{maxrelMCSE}^2 - \mathrm{initES}\).
Numeric positive, default 2: number of initial "burn-in" samples, separated by the Expected Sample Size, to be discarded. Note that the Monte Carlo chain typically starts in a high-probability region, so there is no reason to discard many initial samples.
Integer or NULL (default): thin out the Monte Carlo samples by this value. If NULL: let the diagnostics decide the thinning value.
Logical, default TRUE: output the progress to terminal? If FALSE, the progress is outputted to the file 'main.log' in the outputdir directory.
Logical, default TRUE: save plots of the Monte Carlo traces of diagnostic values?
Logical, default FALSE: save plots of the Monte Carlo traces of the K parameter?
Logical, default FALSE: save plots of the Monte Carlo traces of the Alpha parameter?
List: hyperparameters of the hyperprior; see values in "Usage".
This function takes as main inputs a set of data and metadata, and computes the full joint probability distribution for new data, including its variability. From this full joint distribution any other distributions of interest can subsequently be computed; see Pr() and related functions. This computation can also be interpreted as an estimation of the full joint frequency distribution of the variates in the whole population, beyond the sample data, together with its uncertainty. The computation allows for the use of datapoints with partially missing variables: imputation is automatically made. This imputation is principled, made according to the rules of probability theory.
The output is a "learnt" object, typically saved in a learnt.rds file, which is used in all subsequent probabilistic computations. Other information about the computation is provided in logs and plots, saved in a directory specified by the user.
See vignette('intro') for introductory examples.
The computation is "non-parametric": probability or frequency distributions are not assumed to be Gaussian or of any other specific shape; no "model" is assumed. The mathematical representation of the space of joint frequency distributions follows ideas of Dunson & Bhattacharya (2011); see technical manual for details.
The computation is done via Markov-chain Monte Carlo, using the package Nimble. "Convergence" of the Monte Carlo computation is automatically assessed with methods described in Vehtari & al. (2021) and Kwon & al. (2025); see technical manual for details. The default values for convergence require that all of the following three conditions be fulilled:
The computation's numerical error (Monte-Carlo Standard Error) for the posterior probability must be smaller than 4.7% of the standard deviation of the posterior's variability.
The computation's numerical error for the 0.055- and 0.945-percentiles of the posterior's variability should be smaller than 4.7% of the distance between them.
Typically this requirement leads to final results obtained with the Pr() function having at least two significant digits.
The learn() function can take hours or even days to perform its computations, depending on the size of the dataset, number of variates, and the (initially unknown) "shape" of the underlying probability distribution. For this reason it is typically called within an R script, executed via utils::Rscript. For example, a script 'myscript.R' could have the following structure:
library('prova')learn(
data = 'filename_with_data.csv', # CSV file containing the dataset
metadata = 'filename_with_metadata.csv', # CSV file containing the metadata
outputdir = 'some_directory', # path to output directory
parallel = 8 # machine has more than 8 cores, so we use 8
## possibly other arguments to learn()
)
and then be called on a bash terminal with
$ Rscript myscript.R > learnoutput.log 2>&1 &
with such a call, the file 'learnoutput.log' will contain information about how the computation is proceeding and the estimated end time.
For the mathematical representation of the frequency space:
Dunson, Bhattacharya (2011): Nonparametric Bayes regression and classification through mixtures of product kernels doi:10.1093/acprof:oso/9780199694587.003.0005.
Ishwaran, Zarepour (2002): Exact and approximate sum representations for the Dirichlet process doi:10.2307/3315951.
Porta Mana https://github.com/pglpm/prova/raw/main/development/manual/pglpm2024-bayes_nonparam.pdf.
About Bayesian inference under exchangeability ("population inference"):
Lindley, Novick (1981): The role of exchangeability in inference, doi:10.1214/aos/1176345331.
Bernardo, Smith (2000): Bayesian Theory. Wiley doi:10.1002/9780470316870.
Porta Mana https://github.com/pglpm/prova/raw/main/development/manual/pglpm2024-bayes_nonparam.pdf.
About nonparametrics:
Müller et al. (2015): Nonparametric Bayesian inference. IMS doi:10.1007/978-3-319-18968-0.
Hjort et al. (2010): Bayesian Nonparametrics. Cambridge University Press doi:10.1017/CBO9780511802478.
About Markov-chain Monte Carlo and "convergence":
de Valpine, Paciorek, Turek, & al. (2026): NIMBLE: MCMC, Particle Filtering, and Programmable Hierarchical Modeling doi:10.5281/zenodo.1211190, https://cran.r-project.org/package=nimble.
Kwon & al. (2025): MCMC stopping rules in latent variable modelling doi:10.1111/bmsp.12357.
Vehtari & al. (2021): Rank-normalization, folding, and localization: an improved R-hat for assessing convergence of MCMC doi:10.1214/20-BA1221.
Roy (2020): Convergence diagnostics for Markov chain Monte Carlo doi:10.1146/annurev-statistics-031219-041300.
Gilks & al. (1998): Markov Chain Monte Carlo in Practice. Chapman & Hall/CRC doi:10.1201/b14835.
D. J. C. MacKay (2005): Information Theory, Inference, and Learning Algorithms. Cambridge University Press https://www.inference.org.uk/itila/book.html.
Porta Mana https://github.com/pglpm/prova/raw/main/development/manual/pglpm2024-bayes_nonparam.pdf.
metadatatemplate() to help writing metadata files.
Pr() to calculate probabilities, and qPr() to calculate quantiles, given the data processed by learn().
rPr() to generate datapoints similar to the data processed by learn().
mutualinfo() to calculate mutual information given the data processed by learn().
pread.csv() and pwrite.csv() to read and write CSV files in the format used by learn().
# \donttest{
### WARNING: the following example, if run, might even take a minute or more.
## Create dataset with 3 points of variate 'V' for demonstration:
dataset <- data.frame(V = rnorm(n = 3))
## Create metadata file:
metadata <- data.frame(name = 'V', type = 'continuous')
## Learn from the data:
learnt <- learn(
data = dataset, metadata = metadata,
## the following parameters are unrealistic
## only used to reduce computation time for this example
nsamples = 10, nchains = 1,
startupMCiterations = 10, maxMCiterations = 10,
minESS = 0, initES = 0
)
## Check structure of `learnt` object:
str(learnt)
# }
Run the code above in your browser using DataLab