DREAM Sampling Templates

Collaboration diagram for DREAM Sampling Templates:

Typedefs
using	TasDREAM::DreamDomain = std::function< bool(std::vector< double > const &x)>
	Generic test function whether a sample belongs in the domain. More...
using	TasDREAM::DreamPDF = std::function< void(const std::vector< double > &candidates, std::vector< double > &values)>
	Generic probability distribution used by Tasmanian. More...
using	TasDREAM::DreamModel = std::function< void(const std::vector< double > &candidates, std::vector< double > &outputs)>
	Generic model signature used by Tasmanian. More...
using	TasDREAM::DreamLikelihood = std::function< void(TypeSamplingForm form, const std::vector< double > &model_outputs, std::vector< double > &likely)>
	Generic likelihood signature used by Tasmanian. More...
using	TasDREAM::DreamPrior = std::function< void(TypeSamplingForm form, const std::vector< double > &candidates, std::vector< double > &values)>
	Generic signature for the prior distributions used by Tasmanian. More...
using	TasDREAM::DreamMergedLikelyModel = std::function< void(const std::vector< double > &candidates, std::vector< double > &values)>
	Generic signature for a combination of a likelihood and a model. More...

Functions
DreamDomain	TasDREAM::hypercube (std::vector< double > const &lower, std::vector< double > const &upper)
	Make a lambda that matches the inside signature in SampleDREAM(), test if the vector x is in the hyperbube described by lower and upper.
void	TasDREAM::no_update (std::vector< double > &)
	Dummy function that does not make any changes to the vector as default for the independent_update() in SampleDREAM(). More...
double	TasDREAM::const_one ()
	Dummy function that returns 1.0, used as default for the differential_update() in SampleDREAM(). More...
template<int weight_percent>
double	TasDREAM::const_percent ()
	Template that returns a constant based on the percentage, i.e., weight_percent / 100.0. More...
void	TasDREAM::uniform_prior (TypeSamplingForm, const std::vector< double > &, std::vector< double > &values)
	Uniform prior distribution for both regular and log form. More...
template<TypeSamplingForm form = regform>
DreamPDF	TasDREAM::posterior (DreamModel model, DreamLikelihood likelihood, DreamPrior prior)
	Combines the three components of a Bayesian posterior into a single distribution. More...
template<TypeSamplingForm form = regform>
DreamPDF	TasDREAM::posterior (DreamMergedLikelyModel likelihood_model, DreamPrior prior)
	Overload where the model and likelihood are combined into a single call. More...
template<TypeSamplingForm form = regform>
void	TasDREAM::SampleDREAM (int num_burnup, int num_collect, DreamPDF probability_distribution, DreamDomain inside, TasmanianDREAM &state, std::function< void(std::vector< double > &x)> independent_update=no_update, std::function< double(void)> differential_update=const_one, std::function< double(void)> get_random01=tsgCoreUniform01)
	Core template for the sampling algorithm. More...
template<TypeSamplingForm form = regform>
void	TasDREAM::SampleDREAM (int num_burnup, int num_collect, DreamPDF probability_distribution, DreamDomain inside, TasmanianDREAM &state, TypeDistribution dist, double magnitude, std::function< double(void)> differential_update=const_one, std::function< double(void)> get_random01=tsgCoreUniform01)
	Overload of SampleDREAM() assuming independent update from a list of internally implemented options. More...

Detailed Description

Templates and auxiliary methods for the DREAM sampling. The main template is TasDREAM::SampleDREAM() with one overload and several helper functions. The helpers provide ways to define the probability distribution: either custom defined, interpolated with a sparse grid, or product of Bayesian inference problem.

Typedef Documentation

DreamDomain

using TasDREAM::DreamDomain = typedef std::function<bool(std::vector<double> const &x)>

Generic test function whether a sample belongs in the domain.

The function accepts a single vector with size equal to the number of domain dimensions and must return true if the vector belongs to the domain or false otherwise.

DreamPDF

using TasDREAM::DreamPDF = typedef std::function<void(const std::vector<double> &candidates, std::vector<double> &values)>

Generic probability distribution used by Tasmanian.

The probability distribution must be set to accept multiple candidates and return the value of the unscaled probability distribution for each point.

Parameters

candidates	is a vector with size that is an even multiple of the dimensions, the vector is organized logically into strips of size num-dimensions. Each strip represents a test point that is guaranteed to be inside the domain.
values	is a vector with size equal to the number of strips (samples) in candidates. The vector should not be resized, instead each value has to be overwritten with the corresponding unscaled probability distribution.

Note: the generic probability distribution does not accept a parameter to specify the sampling form, e.g., TasDREAM::logform. The returned values must corresponding to the sampling form set in the TasDREAM::SampleDREAM() template. If logarithm form is used, the values can be negative, in regular form the values must be positive, Tasmanian will not throw an exception but using negative values with TasDREAM::regform leads to undefined behavior.

DreamModel

using TasDREAM::DreamModel = typedef std::function<void(const std::vector<double> &candidates, std::vector<double> &outputs)>

Generic model signature used by Tasmanian.

The model is very similar to the TasDREAM::DreamPDF, in fact the input candidates is the same. The differences are two:

• the model may have multiple outputs and the number of returned outputs must match the number of outputs used by the likelihood.
• the outputs vector will not be set to the correct size, it must be resized, because Tasmanian does not keep track of the number of outputs.

Note: the outputs vector will be fed as input to the TasDREAM::DreamLikelihood.

DreamLikelihood

using TasDREAM::DreamLikelihood = typedef std::function<void(TypeSamplingForm form, const std::vector<double> &model_outputs, std::vector<double> &likely)>

Generic likelihood signature used by Tasmanian.

The likelihood assigns a value to how likely the model outputs are given data (i.e., measurements). Classes that inherit TasmanianLikelihood will automatically convert to a lambda object with this signature.

Parameters

form	is the sampling type used in the call to TasDREAM::SampleDREAM(), in a custom object it is sufficient to implement only one form, Tasmanian likelihood classes implement both and thus the parameter specifies which implementation to use.
model_outputs	is identical to outputs parameter in TasDREAM::DreamModel except here it is used as an input to compute the likelihood.
likely	is a vector with size equal to the entries with outputs, the entries must be overwritten with the corresponding likelihood.

DreamPrior

using TasDREAM::DreamPrior = typedef std::function<void(TypeSamplingForm form, const std::vector<double> &candidates, std::vector<double> &values)>

Generic signature for the prior distributions used by Tasmanian.

Specifies the prior distribution, the TasDREAM::uniform_prior() satisfies this signature.

Parameters

form	is the same as in TasDREAM::DreamLikelihood.
candidates	is the same as in TasDREAM::DreamModel and TasDREAM::DreamPDF.
values	are similar to the likely in TasDREAM::DreamLikelihood, but instead of the likelihood the values define the prior distribution in the specified form.

DreamMergedLikelyModel

using TasDREAM::DreamMergedLikelyModel = typedef std::function<void(const std::vector<double> &candidates, std::vector<double> &values)>

Generic signature for a combination of a likelihood and a model.

The likelihood and the model are not always separated, e.g., a sparse grid approximation can be used to interpolated the likelihood which has a single output and is therefore cheaper than interpolating multi-output model. The implementation should be equivalent to:

Parameters

candidates	is same as in TasDREAM::DreamModel.
values	is same as in TasDREAM::DreamLikelihood likely parameter.

Function Documentation

no_update()

void TasDREAM::no_update ( std::vector< double > &  )

inline

Dummy function that does not make any changes to the vector as default for the independent_update() in SampleDREAM().

The function is no-op.

const_one()

double TasDREAM::const_one ( )

inline

Dummy function that returns 1.0, used as default for the differential_update() in SampleDREAM().

Just an inline function that returns 1.0.

const_percent()

template<int weight_percent>

double TasDREAM::const_percent

(

)

Template that returns a constant based on the percentage, i.e., weight_percent / 100.0.

The template simplifies the syntax when calling SampleDREAM() with a constant differential update. For example, setting the update to 0.5 can be done with

TasmanianDREAM state(...);
SampleDREAM(..., independent_update, const_percent<50>, state);

uniform_prior()

void TasDREAM::uniform_prior ( TypeSamplingForm  , const std::vector< double > &  , std::vector< double > &  values  )

inline

Uniform prior distribution for both regular and log form.

Applies uniform (non-informative) prior, can be used with any of the Bayesian inference methods. In practice, this actually does nothing, since adding zero (in logform) or mulitplying by 1 (in regform) amounts to nothing.

posterior() [1/2]

template<TypeSamplingForm form = regform>

DreamPDF TasDREAM::posterior	(	DreamModel	model,
		DreamLikelihood	likelihood,
		DreamPrior	prior
	)

Combines the three components of a Bayesian posterior into a single distribution.

The Bayesian posterior has a model, likelihood function and a prior distribution. This function combines three function objects into a single probability distribution that be passed to SampleDREAM, e.g.,

SampleDREAM(num_burnup, num_collect, posterior(likely, model, uniform_prior), ...);

The generality of the approach used here comes at the price of volatility. There is no builtin error-checking and error-detection on vector sizes. Specifically, the number of inputs provided by the model must match the outputs accepted by the likelihood, and the number of dimensions accepted by the model and prior must be the same.

Template Parameters

form

indicates whether to use the regular or logarithm of the sampling problem; Gaussian-types of likelihood functions are often used where the regular form can be while the log-form is , working with a simple quadratic function can be more stable with respect to rounding error. The template parameter must be the same as in the call to TasDREAM::SampleDREAM().

Parameters

model	accepts a set of model inputs and will return the corresponding model values. candidates is the same as in the input of probability_distribution() in the call to TasDREAM::SampleDREAM(). Logically the candidates will be arranged in strips of size equal to the problem dimensions, the vector size will divide evenly by the dimensions and the factor is the number of candidates. outputs must be resized to match the number of candidates times the number of outputs, the behavior must match that of TasmanianSparseGrid::evaluateBatch().
likelihood	accepts a set of model outputs and provides a measure of how likely those outputs are given some observed data with some noise. Tasmanian provides likelihood functions that can be used here, e.g., TasDREAM::LikelihoodGaussIsotropic and TasDREAM::LikelihoodGaussAnisotropic. The TasDREAM::TypeSamplingForm will always match the template parameter form, thus, it is sufficient to implement only one sampling from. The Tasmanian likelihood classes implement both forms, hence the added flexibility. The model_outputs is a vector with size equal to the number of candidates times the number of outputs, i.e., must match the output of the model. The likely will have size equal to the number of candidates and must be filled (without resize) with the likelihood for each set of model outputs.
prior	provides the values of the prior distribution in either regular or logarithm form. The prior will take in the same candidates as the model and a vector of the same size as likely, and must return the values of the corresponding prior distribution in either regular or logarithm form.

Example usage:

auto model = TasGrid::read("foo");
TasDREAM::LikelihoodGaussIsotropic likely(0.1, data);
TasDREAM::SampleDREAM(..., TasDREAM::posterior(likely, model, TasDREAM::uniform_prior), ...);

posterior() [2/2]

template<TypeSamplingForm form = regform>

DreamPDF TasDREAM::posterior	(	DreamMergedLikelyModel	likelihood_model,
		DreamPrior	prior
	)

Overload where the model and likelihood are combined into a single call.

There are situations where splitting the model and likelihood is undesirable, e.g., if the model has a huge number of outputs it may be easier to construct a sparse grid surrogate to the single-output combined model and likelihood. This is a short hand-template that uses such model-likelihood and combines it with a prior distribution.

SampleDREAM() [1/2]

template<TypeSamplingForm form = regform>

void TasDREAM::SampleDREAM	(	int	num_burnup,
		int	num_collect,
		DreamPDF	probability_distribution,
		DreamDomain	inside,
		TasmanianDREAM &	state,
		std::function< void(std::vector< double > &x)>	independent_update = no_update,
		std::function< double(void)>	differential_update = const_one,
		std::function< double(void)>	get_random01 = tsgCoreUniform01
	)

Core template for the sampling algorithm.

Evolves the chains of the state using the Metropolis algorithm where the updates are comprised of two components, independent and differential. The independent component relies on independently sampled (pesudo)-random numbers with zero mean (could be deterministic constant zero). The differential component is based on the difference between two randomly chosen chains, which effectively exchanges information between the chains.

The implementation is very generic using lambda objects to describe most aspects of the problem, including the probability distribution, domain, etc. However, the generality comes with some sacrifice in resilience, i.e.,

• each lambda object must respect the problem dimensions, the template will populate the inputs with the correct number of entries, but the lambdas have to properly utilize the entries.
• each lambda will have to capture external objects and variables and should avoid capture by value for large vectors, e.g., a sparse grid or a large data vector, while objects captured by reference must remain alive during the call to the template. See couple of examples after the variable listing.

Template Parameters

form

indicates whether the probability_distribution() function return the regular form or the logarithm of the desired pdf; in some cases, taking the exponential of a very large negative number can lead to problems with rounding comparison between numbers very close to zero, hence the logarithm from could be more stable from the round-off error perspective.

Parameters

num_burnup	is the number of initial iterations that will not be saved in the history; the Metropolis algorithm guarantees convergence to the desired distribution but only in the limit; when working with a finite number of samples the results can be contaminated by the initial state which can be significantly different from the limit. A common practice is to take num_burnup to be roughly equal to num_collect, but that is only a heuristic suggestion.
num_collect	is the number of iterations that will be saved in the state history, the total number of collected samples is num_collect times state.getNumChains().
probability_distribution	must accept a vector of candidates locations and return the values of the unscaled probability distribution at those points. The input-output relation is similar to working with TasmanianSparseGrid::evaluateBatch() when the grid has a a single output. The candidates input will be logically divided into strips of size state.getNumDimensions(), the total number of strips will be between 1 and state.getNumChains() depending on the number of proposals that fall within the domain, i.e., the inputs will always satisfy the inside() condition. The values will have size equal to the number of strips and each entry should be filled the corresponding value of the probability density function. The values vector should not be resized. Any TasGrid::TasmanianSparseGrid object with a single output and non-zero loaded points can be passed in as a probability_distribution; when using TasDREAM::regform the output of the grid must be always non-negative, using TasDREAM::logform has no lower bound requirements. In the context of Bayesian inference, there probability distribution is comprised of three components: model, likelihood, and prior. See TasDREAM::posterior().
inside	must accept a vector of size equal tot he number of dimensions and return true if the vector describes a point in the domain and false otherwise. The function will be called for every proposed sample and unlike the probability_distribution only one vector will be given as input at a time. Each TasGrid::TasmanianSparseGrid object some with a canonical domain that depends on the rule and can be optionally transformed, the TasmanianSparseGrid::getDomainInside() method will produce a lambda object describing the sparse grid domain. See also TasDREAM::hyperbube().
state	must be an initialized instance of TasmanianDREAM. The number of chains and dimensions must be set to match the lambda objects and TasmanianDREAM::setState() must be called to load the state with a valid initial set of vectors. The state will be evolved in total of num_burnup + num_collect iterations and the last set of num_collect iterations will be appended to the state history.
independent_update	must accept a vector of size state.getNumDimensions() and (without resizing) perturb each entry by adding an independent zero-mean random number. It is possible to omit the independent_update, e.g., pass an object that will not modify the input vector x; however, this can lock the chains to a set of values not dense in the domain which in turn will break convergence. For example, if the domain is the real line, the initial state has integer entries, and the differential_update is set to 100%, then all proposals will be integers, non-integer values will never be considered. TasDREAM::SampleDREAM provides an overload where the independent update is set to a distribution from an included list with known magnitude.
differential_update	is a random or deterministic number between zero and one, indicating the magnitude of the difference between two randomly chosen chains that will be added to the next proposal. Using deterministic zero will decouple the chains, i.e., the state of one chain will never affect any other chain and the algorithm will reduce to the propagation of multiple independent chains of the standard Metropolis-Hastings method with the selected independent_update. Using negative values (by symmetry) is equivalent to using a positive value with the same magnitude, and values larger than 1 are likely to result in poor mixing and bad convergence rate. The default differential update is deterministic one and deterministic percentage can be specified with the const_percent() template.
get_random01	is the pseudo-random number generator to be used in the sampling procedure. By default, Tasmanian will use rand() divided by RAND_MAX, but this is implementation dependent and not always optimal.

Correct call using a sparse grid object as input:

auto grid = TasGrid::read("foo"); // create a grid object
SampleDREAM(..., grid, ...);
// above, the grid will create a lambda object that will capture grid by reference
// the lambda object is destroyed at the end of the call and before grid

Incorrect call, never do that:

SampleDREAM(..., TasGrid::read("foo"), ...); // <- Incorrect
// above, read() will create a TasmanianSparseGrid object which will create a lambda
// but the grid will be destroyed before entering the SampleDREAM() and cause a segfault.

SampleDREAM() [2/2]

template<TypeSamplingForm form = regform>

void TasDREAM::SampleDREAM	(	int	num_burnup,
		int	num_collect,
		DreamPDF	probability_distribution,
		DreamDomain	inside,
		TasmanianDREAM &	state,
		TypeDistribution	dist,
		double	magnitude,
		std::function< double(void)>	differential_update = const_one,
		std::function< double(void)>	get_random01 = tsgCoreUniform01
	)

Overload of SampleDREAM() assuming independent update from a list of internally implemented options.

Uses independent update is applied uniformly to all dimensions and comes from a list of internal functions, e.g., uniform or Gaussian. This overload wraps around functions such as TasDREAM::applyUniformCorrection() and TasDREAM::applyGaussianCorrection().