Markov chain Monte Carlo Search Algorithm¶
markov_chain_monte_carlo¶
Module containing algorithm classes for Markov chain Monte Carlo sampling.
The McMC algorithms are initialised using the McMCAlgorithmCreator class:
-
class
MTfit.algorithms.markov_chain_monte_carlo.McMCAlgorithmCreator[source]¶ Generates the Markov chain Monte Carlo algorithm from parameters
Creates the correct McMC algorithm object depending on the selected mode.
McMC Creator
Returns initialised object of desired type. This is extensible using the MTfit.algorithms.markov_chain_monte_carlo group in pkg_resources (see extensions documentation) to add additional options
- Args
mode:[MetropolisHastings] McMC mode to use options are: MetropolisHastings alpha:[False] Default parameters for Markov chain steps, depend on transition probabilities trans_dimensional:[False] Boolean flag to run with trans-dimensional sampling learning_length:[10000] Number of samples for learning period and to discard from chain. parameterisation:['Tape'] Source type parameterisation to use options are: Tape transition:['Gaussian'] Transition pdf to use options are: 'Gaussian' chain_length:[1000000] End point of the Markov chain. min_acceptance_rate:[0.3] Minimum targetted sample acceptance rate. max_acceptance_rate:[0.5] Maximum targetted sample acceptance rate. acceptance_rate_window:[1000] Number of samples to use in learning period for calculating and
modifying the acceptance rate.initial_sample:['Grid'] Initialisation sampling mode to use options are 'Grid'.
- Returns
- McMC sampling object.
-
static
__new__(self, mode='metropolis_hastings', alpha=False, trans_dimensional=False, learning_length=10000, parameterisation='tape', transition='gaussian', chain_length=1000000, min_acceptance_rate=0.3, max_acceptance_rate=0.5, acceptance_rate_window=1000, initial_sample='grid', **kwargs)[source]¶ McMC Creator
Returns initialised object of desired type. This is extensible using the MTfit.algorithms.markov_chain_monte_carlo group in pkg_resources (see extensions documentation) to add additional options
- Args
mode:[MetropolisHastings] McMC mode to use options are: MetropolisHastings alpha:[False] Default parameters for Markov chain steps, depend on transition probabilities trans_dimensional:[False] Boolean flag to run with trans-dimensional sampling learning_length:[10000] Number of samples for learning period and to discard from chain. parameterisation:['Tape'] Source type parameterisation to use options are: Tape transition:['Gaussian'] Transition pdf to use options are: 'Gaussian' chain_length:[1000000] End point of the Markov chain. min_acceptance_rate:[0.3] Minimum targetted sample acceptance rate. max_acceptance_rate:[0.5] Maximum targetted sample acceptance rate. acceptance_rate_window:[1000] Number of samples to use in learning period for calculating and
modifying the acceptance rate.initial_sample:['Grid'] Initialisation sampling mode to use options are 'Grid'.
- Returns
- McMC sampling object.
The standard McMC algorithms used are the multiple try algorithms using the Tape and Tape (2012) parameterisation and a Gaussian proposal function.
-
class
MTfit.algorithms.markov_chain_monte_carlo.IterativeMultipleTryMetropolisHastingsGaussianTape(*args, **kwargs)[source]¶ Metropolis-Hastings Markov chain Monte Carlo Algorithm using Gaussian Transition PDF and Tape parameterisation, chain length ends maximum length
Algorithm ends when the number of samples in the chain equals the chain_length, where the chain length corresponds to the number of accepted samples. The parameterisation is from Tape and Tape (A geometric setting for moment tensors, Tape and Tape, 2012, GJI 190 pp 476-490).
This is a child class of the IterativeMetropolisHastingsGaussianTape
Initialisation of Metropolis-Hastings Markov chain Monte Carlo Algorithm using Gaussian Transition PDF and Tape parameterisation, constrained by chain length
The chain length corresponds to the number of accepted samples.
- Keyword Args
- number_samples:[1000] Maximum number of samples to try on each iteration. (unused samples are dropped)
This is a child class of the IterativeMetropolisHastingsGaussianTape
-
acceptance(x, ln_likelihood_x, *args)¶ Calculates acceptance
Calculates the acceptance from the Metropolis condition.
- Args
- x: Model values ln_likelihood_x: Model ln_likelihood.
- Returns
- float:acceptance
-
acceptance_rate()¶ Gets acceptance rate.
- Returns
- float: acceptance rate.
-
convert_sample(x)¶ Converts sample to MT form
- Args
- x: sample
- Returns
- Converted Sample
-
eigenvectors_mt_2_mt6(diag, a, b, c)¶ Converts eigenvectors and eigenvalues to moment tensor 6-vector.
Generates moment tensor 6-vector.
- Args
- diag: numpy matrix of eigenvalues. a: numpy array of eigenvectors corresponding to the first eigenvalue. b: numpy array of eigenvectors corresponding to the second eigenvalue. c: numpy array of eigenvectors corresponding to the third eigenvalue.
- Returns
- numpy matrix of moment tensor 6-vectors.
-
get_sampling_model(kwargs, file_sample, file_safe)¶ Get the sampling model from the entry points
-
initialise()¶ Initialse samples
Initialises the chain either using an initialiser if set, otherwise using random_sample.
- Returns
- new_sample,False
-
is_dc(x)¶ Checks if a sample x is double-couple
-
iterate(result)¶ Iterate from result
- Args
- result: Result dictionary from forward task (e.g. MTfit.inversion.ForwardTask)
- Returns
- new_sample,End where End is a boolean flag to end the chain if the length of accepted samples is longer than the chain length.
-
learning_check()¶ Check if still in learning period
-
new_sample(jump=0.0, gaussian_jump=False)[source]¶ Get new samples including multiple samples to try
-
output(*args, **kwargs)¶ Returns output dictionary
- Returns
- dict: Output dictionary
-
prior(x, dc=None, basic_cdc=None, *args, **kwargs)¶ Evaluates prior probability for x
- Returns
- float: prior probability
-
random_basic_cdc()¶ Generates a random C+DC
Generates a random C+DC using flat pdfs for each parameter (gamma,delta given by alpha).
- Returns
- Random MT in Tape parameterisation
-
random_clvd(*args)¶ Generate random CLVD moment tensors (size 6,number_samples)
Generates random CLVD moment tensors with random orientations.
- Returns
- numpy matrix of random CLVD moment tensors, size 6,number_samples
-
random_dc()¶ Generates a random DC
Generates a random DC using flat pdfs for each parameter (gamma=delta=0).
- Returns
- Random MT in Tape parameterisation
-
random_orthogonal_eigenvectors()¶ Generates random orthogonal eigenvectors.
- Returns
- list of numpy arrays of eigenvectors.
-
random_sample()¶ Return random samples
Returns: results from generating the random samples.
-
random_type(diag)¶ Generate random moment tensors from a given set of eigenvalues.
Generates random moment tensors with a given set of eigenvalues and random orientations.
- Args
- diag: numpy matrix of eigenvalues.
- Returns
- numpy matrix of random moment tensors, size 6,number_samples
-
transition_pdf(x, x1, dc=None, basic_cdc=None)¶ Calculates gaussian transition pdf
Evaluates the gaussian transition probability for x and x1
- Returns
- float: transition probability
-
class
MTfit.algorithms.markov_chain_monte_carlo.IterativeMultipleTryTransDMetropolisHastingsGaussianTape(*args, **kwargs)[source]¶ Trans-Dimensional Metropolis-Hastings Markov chain Monte Carlo Algorithm using Gaussian Transition PDF and Tape parameterisation, chain length ends maximum length
Trans Dimensional sampling, jumping between double-couple and full moment tensor models. Algorithm ends when the number of samples in the chain equals the chain_length, where the chain length corresponds to the number of accepted samples. The parameterisation is from Tape and Tape (A geometric setting for moment tensors, Tape and Tape, 2012, GJI 190 pp 476-490).
This is a child class of both the IterativeTransDMetropolisHastingsGaussianTape and IterativeMultipleTryMetropolisHastingsGaussianTape
Initialisation of Trans-Dimensional Metropolis-Hastings Markov chain Monte Carlo Algorithm using Gaussian Transition PDF and Tape parameterisation, chain length ends maximum length
This is a child class of both the IterativeTransDMetropolisHastingsGaussianTape and IterativeMultipleTryMetropolisHastingsGaussianTape
-
acceptance(x, ln_likelihoodx, dc_prior=0.5)¶ Calculates acceptance
Calculates the acceptance from the Trans-Dimensional Metropolis condition.
- Args
- x: Model values ln_likelihoodx: Model likelihood.
- Returns
- float:acceptance
-
acceptance_rate()¶ Gets acceptance rate.
- Returns
- float: acceptance rate.
-
convert_sample(x)¶ Converts sample to MT form
- Args
- x: sample
- Returns
- Converted Sample
-
eigenvectors_mt_2_mt6(diag, a, b, c)¶ Converts eigenvectors and eigenvalues to moment tensor 6-vector.
Generates moment tensor 6-vector.
- Args
- diag: numpy matrix of eigenvalues. a: numpy array of eigenvectors corresponding to the first eigenvalue. b: numpy array of eigenvectors corresponding to the second eigenvalue. c: numpy array of eigenvectors corresponding to the third eigenvalue.
- Returns
- numpy matrix of moment tensor 6-vectors.
-
get_sampling_model(kwargs, file_sample, file_safe)¶ Get the sampling model from the entry points
-
initialise()¶ Initialse samples
Initialises the chain either using an initialiser if set, otherwise using random_sample.
- Returns
- new_sample,False
-
is_dc(x)¶ Checks if a sample x is double-couple
-
iterate(*args, **kwargs)¶ Iterates new sample
Calculates the acceptance and handles the pDC output
-
jump_params(x=False)¶ Calculates dimension jump parameters
Calculates either the probabilities of the balancing parameters if x exists, else returns the two balancing parameters.
- Args
- x:[False] Can be given as the two balancing parameters
- Returns
- if x is not False:
- float: q - probability of balancing parameters
- else:
- float: gamma - Randomly distributed gamma value float: delta - Randomly distributed delta value
-
learning_check()¶ Check if still in learning period
-
new_sample()[source]¶ Generates new sample
Generates a new sample by drawing a new sample from the gaussian transition pdf, conditional on the previous sample. Probability if the sample jumping between models given by the dimension_jump_prob kwarg on initialisation.
- Returns
- New sample.
-
output(*args, **kwargs)¶ Return output dict
Return output dict including pDC value
-
prior(x, dc=None, basic_cdc=None, *args, **kwargs)¶ Evaluates prior probability for x
- Returns
- float: prior probability
-
random_basic_cdc()¶ Generates a random C+DC
Generates a random C+DC using flat pdfs for each parameter (gamma,delta given by alpha).
- Returns
- Random MT in Tape parameterisation
-
random_clvd(*args)¶ Generate random CLVD moment tensors (size 6,number_samples)
Generates random CLVD moment tensors with random orientations.
- Returns
- numpy matrix of random CLVD moment tensors, size 6,number_samples
-
random_dc()¶ Generates a random DC
Generates a random DC using flat pdfs for each parameter (gamma=delta=0).
- Returns
- Random MT in Tape parameterisation
-
random_orthogonal_eigenvectors()¶ Generates random orthogonal eigenvectors.
- Returns
- list of numpy arrays of eigenvectors.
-
random_sample()¶ Return random samples
Returns: results from generating the random samples.
-
random_type(diag)¶ Generate random moment tensors from a given set of eigenvalues.
Generates random moment tensors with a given set of eigenvalues and random orientations.
- Args
- diag: numpy matrix of eigenvalues.
- Returns
- numpy matrix of random moment tensors, size 6,number_samples
-
transition_pdf(x, x1, dc=None, basic_cdc=None)¶ Calculates gaussian transition pdf
Evaluates the gaussian transition probability for x and x1
- Returns
- float: transition probability
-
The McMC sampling algorithms inherit from the BaseAlgorithm class:
-
class
MTfit.algorithms.base.BaseAlgorithm(number_samples=10000, dc=False, quality_check=False, file_sample=False, file_safe=True, generate=False, *args, **kwargs)[source]¶ Base class for algorithms including the base methods required for forward modelling based inversions.
BaseAlgorithm initialisation
- Args
number_samples:[10000] Number of samples to use per iteration. dc:[False] Boolean to select inversion constrained to double-couple
space or over the full moment tensor space.- quality_check:[False] Carry out quality check after 40,000 samples.
- If there are too many non-zero samples, the inversion is stopped for that event.
- file_sample:[False] Boolean to select whether to use file sampling
- (saving the output to file).
- file_safe:[False] Boolean to select whether to use file safe outputs
- when using file sampling.
- generate:[False] Boolean to select whether to generate the random
- moment tensors in the forward task or not.
- Keyword Arguments
- basic_cdc:[False] Boolean to select whether to use the Basic CDC
- model (inversion constrained to crack+double-couple space).
- number_events:[1] Integer giving the number of events in this
- inversion (for use with multiple event inversions).
- sampling_prior:[6sphere_prior] string selector for the full moment
- tensor sampling prior (MTfit.sampling_prior entry point)
- sampling:[6sphere] string selector for the full moment tensor
- sampling model (MTfit.sampling entry point)
- model:[False] selector for alternate models using the
- MTfit.sample_distribution entry_point
-
eigenvectors_mt_2_mt6(diag, a, b, c)[source]¶ Converts eigenvectors and eigenvalues to moment tensor 6-vector.
Generates moment tensor 6-vector.
- Args
- diag: numpy matrix of eigenvalues. a: numpy array of eigenvectors corresponding to the first eigenvalue. b: numpy array of eigenvectors corresponding to the second eigenvalue. c: numpy array of eigenvectors corresponding to the third eigenvalue.
- Returns
- numpy matrix of moment tensor 6-vectors.
-
get_sampling_model(kwargs, file_sample, file_safe)[source]¶ Get the sampling model from the entry points
-
initialise()[source]¶ Basic initialisation function
Return the first task
- Returns
- task, False task: New task for forward model
-
iterate(result)[source]¶ Basic iteration function
- Args
- result: Result from forward Task
- Returns
- task,End task: New task for forward model End: Boolean flag for whether the inversion is finished.
-
output(normalise=True, convert=False, discard=10000)[source]¶ Return the algorithm results for output.
- Returns
- Dictionary containing output.
-
random_clvd(*args)[source]¶ Generate random CLVD moment tensors (size 6,number_samples)
Generates random CLVD moment tensors with random orientations.
- Returns
- numpy matrix of random CLVD moment tensors, size 6,number_samples
-
random_dc(*args)[source]¶ Generate random double-couple moment tensors (size 6, number_samples)
Generates random double-couple moment tensors with random orientations.
- Returns
- numpy matrix of random double-couple moment tensors, size 6,number_samples