Sam Class

class sam.Sam(logProbability, scale, lowerBounds=None, upperBounds=None)

A class for sampling from probability distributions.

Instantiates the sampler class and sets the logProbability function.

Parameters:
  • logProbability – A function which takes one or three arguments and returns the natural log of the probability evaluated at the position given by the first argument. If three arguments are allowed, the second argument provides an array to write the gradient to, and the third argument is a bool indicating whether a gradient is required for that call to the function.
  • scale – An array whose length is the number of parameters in the target distribution.
  • lowerBounds – An array whose length is the number of parameters which defines lower boundaries below which the parameters will not be sampled. The sampler is optimized so that these boundaries are enforced efficiently, and so will not decrease the acceptance rate. If None, no boundaries are enforced.
  • upperBounds – Same as lowerBounds, but defines the upper boundaries
Returns:

An instantiated object of the class.

Running the Sampler

sam.Sam.run(self, Size nSamples, x, Size burnIn=0, Size thinning=0, Size recordStart=0, Size recordStop=-1, bool collectStats=False, Size threads=1, bool showProgress=True)

Begin sampling the parameters from the given logProbability dist.

Parameters:
  • nSamples – The desired number of samples to record per thread.
  • x – The initial position(s) to start the sampler at.
  • burnIn – The number of MCMC steps to take before recording begins.
  • thinning – The number of MCMC steps to take between recordings and burn-in steps. This directly multiplies the amount of work the sampler needs to do.
  • recordStart – The index of the first parameter to be recorded. Default is 0.
  • recordStop – The index of the last parameter to be recorded, plus one. Default is the number of parameters (all are recorded).
  • collectStats – Whether the sampler should collect running statistics as it runs. This is probably only desirable if not all parameters are being recorded. Default is False
  • threads – The number of computational threads to run in. If this is greater than 1, the multiprocessing library is used to run this many fully independent samplers in parallel.
  • showProgress – Whether to print a progress bar as the sampler runs.
Returns:

The parameter samples, of the shape [N x M X T], where N is the number of parameters of target distribution, M is the number of samples requested (nSamples), and T is the number of threads used. If threads is 1 the last dimension is omitted. This data is also stored internally and can be accessed later using other functions.

Picking Samplers

By default, Sam will use a simple Metropolis-Hastings sampler. You can manage the different samplers by using the following functions. Note that you can several samplers in combination, so as to use different algorithms for different parameters. This is done using the dStart and dStop

sam.Sam.addMetropolis(self, covariance=None, Size dStart=0, Size dStop=-1)

Adds a metropolis sampler with a non-diagonal covariance.

This sampler sets up a Metropolis-Hastings sampler to be used during the sampling procedure. The proposal distribution is a normal distribution centered at self.x, with a covariance supplied by the user. If no covariance is supplied, default is a matrix with self.scale (set during initialization) as the diagonal (this is slightly faster to produce random numbers for).

Parameters:
  • covariance – The covariance matrix to be used. Should be [M x M], where M=dStop-dStart. May be None to use a diagonal matrix with self.scale as the diagonals.
  • dStart – The index of the first parameter to be included. Default is zero.
  • dStop – The index of the last parameter to be included, plus one. Default is the last index + 1.
sam.Sam.addAdaptiveMetropolis(self, covariance=None, int adaptAfter=-1, int recordAfter=-1, int refreshPeriod=100, double scaling=-1, double eps=1e-9, Size dStart=0, Size dStop=-1)

Adds an Adaptive Metropolis sampler to the sampling procedure.

This sampler is the Adaptive Metropolis (AM) algorithm presented in Haario et al. (2001). The algorithm initially samples with a given proposal covariance, but after a number of steps, uses the covariance of the samples to estimate the optimal proposal covariance. Each time the propsal is updated, the cholesky (order n^3) must be computed, so it may be advisable not to set the refresh period too low if n is large. Note that the estimated covariance is updated every time the sampler is called, but that this is not propagated to the sampler until the refresh occurs.

Parameters:
  • covariance – The initial proposal covariance to sample with. Should be [M x M], where M = nStop - nStart.
  • adaptAfter – How many times the sampler must be called (and thereby collect samples) before the adapted covariance is used. Must be larger than the number of dimensions being adapted to. The default (triggered by any negative number) is three times that or 100, whichever is greater.
  • recordAfter – How many times the sampler must be called before the adapted covariance begins to take in samples. This is to prevent pre-burned-in samples from being used to adapt with, which can dramatically reduce the effectiveness of sampling.
  • refreshPeriod – How many times the sampler is called between the cholesky of the covariance being updated (the expensive part).
  • scaling – How much to scale the estimated covariance to get the proposal covariance. Default, signified by -1, is to use 5.76/nDim, suggested in Haario et al.
  • eps – The epsilon parameter in Haario et al. (2001). It needs to be small but nonzero for the theory to work, but in practice seems to work at zero as well. Default is 1e-9, and probably will not need to be changed by the user unless the covariance is on that scale or less.
  • dStart – The index of the first parameter to be included. Default is zero.
  • dStop – The index of the last parameter to be included, plus one. Default is the last index + 1.
sam.Sam.addHMC(self, Size nSteps, double stepSize, Size dStart=0, Size dStop=-1)

Adds a Hamiltonian Monte Carlo sampler to the sampling procedure.

This sampler sets up an HMC sampler to be used during the sampling procedure. The proposal distribution is a normal distribution centered at self.x, with diagonal covariance where the scale (set when the sampler object was initialized) is the diagonal components.

Parameters:
  • nSteps – The number of steps in the trajectory
  • stepSize – Multiplied by self.scale to scale the process.
  • dStart – The index of the first parameter to be included. Default is zero.
  • dStop – The index of the last parameter to be included, plus one. Default is the last index + 1.
sam.Sam.printSamplers(self)

Prints the list of any/all sampling systems set up so far.

This refers to samplers added using e.g. addMetropolis() functions. See the overall class documentation for more information on these. The type of function and any parameters given (unique to each sampler type) are printed.

sam.Sam.clearSamplers(self)

Clears the list of samplers.

Surrogate Sampling

Sam supports the use of a Gaussian process (GP) as a surrogate model for the posterior probability. This means that each time the sampler needs to know the value of the posterior probability at a point, the GP is consulted first. If the GP is able to estimate the posterior probability with uncertainty less than a certain tolerance, then this estimate is used and the logProbability function provided by the user is not called. Otherwise logProbability is called and the result is added to the GP for future reference. This approach can reduce the number of likelihood evaluations considerably at the cost of some modest overhead. It becomes less effective for high-dimensional problems. When the surrogate is enabled, you can access it via Sam.surrogate to e.g. retrieve the results of actual likelihood evaluations ( Sam.surrogate.x and Sam.surrogate.y)

sam.Sam.enableSurrogate(self, xInit, yInit, kernel='matern32', tol=1e-2)

Turns on surrogate sampling and initializes the surrogate GP.

In order to get reasonable parameters for the Gaussian process, you must provide at least some initial points to optimize the GP on. More points can help reduce the number of likelihood evaluations, but will slow down GP computation, so don’t overdo it.

Parameters:
  • xInit – The sample positions to initialize the GP with.
  • yInit – The values of the likelihood to initialize the GP with.
  • kernel – The kernel to use in the GP surrogate. Must be differentiable if you want to use the gradient.
  • tol – How much uncertainty in the log likelihood to permit without calling the likelihood function.
sam.Sam.disableSurrogate(self)

Disables the use of a surrogate model if previously enabled.

Examining the Results

sam.Sam.samples = the samples collected as a Numpy ndarray of shape [nThreads, nSamples, nParameters]
sam.Sam.results = same as ``samples``, but flattend across threads: [(nThreads x nSamples), nParameters]
sam.Sam.summary(self, paramIndices=None, returnString=False)

Prints/returns some summary statistics of the previous sampling run.

Statistics are the parameter index, the acceptance rate, mean, and standard deviation, as well as the 16th, 50th, and 84th percentiles.

Parameters:paramIndices – The indices of the parameters to be described. If set to None, then all parameters will be described.
Returns:The summary message as a string if returnString is True, otherwise returns None.
sam.Sam.getStats(self)

Returns running-average and standard deviation statistics.

Note that this is only applicable if collectStats was enabled during the sampling procedure (default is off). You will need to compute the values directly from the sample otherwise. This is intended as a tool for tracking nuisance parameters for which it was not worth recording the samples from – see recordStart and recordStop in the self.run arguments.

Returns:An array of means and an array of standard deviations of the parameter samples.

Model Comparison

Several functions are available for computing model comparison statistics on the results of a run. Since only the posterior probability is known to the sampler, you must provide a function which returns the prior probability given the parameters in order to compute these.

sam.Sam.getAIC(self, prior)

Approximates the AIC from the sampled values.

This function needs to know the prior because only the full posterior probability was given to the sampler at runtime. The value of the prior is subtracted off, and the sample with the maximum resulting likelihood is identified. The BIC is computed assuming that this is the global maximum likelihood; for a well sampled posterior, this will be a good approximation.

This function shouldn’t be used in a heirarchical setting, as the AIC is not defined for that case (the number of parameters isn’t clearly defined). Consider the DIC instead.

Parameters:prior – a function which takes an array of parameters (same length as the logProbability function) and returns the prior for them.
Returns:The estimated AIC for the model as a floating point number.
sam.Sam.getBIC(self, prior, nPoints)

Approximates the BIC from the sampled values.

This function needs to know the prior because only the full posterior probability was given to the sampler at runtime. The value of the prior is subtracted off, and the sample with the maximum resulting likelihood is identified. The BIC is computed assuming that this is the global maximum likelihood; for a well sampled posterior, this will be a good approximation.

This function shouldn’t be used in a heirarchical setting, as the BIC is not defined for that case (the number of parameters isn’t clearly defined). Consider the DIC instead.

Parameters:
  • prior – a function which takes an array of parameters (same length as the logProbability function) and returns the prior for them.
  • nPoints – The number of data points used to evaluate the likelihood.
Returns:

The estimated BIC for the model as a floating point number.
sam.Sam.getDIC(self, prior)

Approximates the DIC from the sampled values.

This function needs to know the prior because only the full posterior probability was given to the sampler at runtime. The value of the prior is subtracted off, and the sample with the maximum resulting likelihood is identified. The BIC is computed assuming that this is the global maximum likelihood; for a well sampled posterior, this will be a good approximation.

Parameters:prior – a function which takes an array of parameters (same length as the logProbability function) and returns the prior for them.
Returns:The estimated DIC for the model as a floating point number.