bayesbay.discretization.Voronoi

class bayesbay.discretization.Voronoi(name, spatial_dimensions, vmin, vmax, perturb_std, n_dimensions=None, n_dimensions_min=2, n_dimensions_max=10, n_dimensions_init_range=0.3, parameters=None, birth_from='neighbour')

Utility class for Voronoi tessellation

Parameters:

• name (str) – name attributed to the Voronoi tessellation, for display and storing purposes

• spatial_dimensions (int) – number of dimensions of the desired Voronoi tessellation, e.g. 1D, 2D, or 3D.

• vmin (Union[Number, np.ndarray]) – minimum/maximum value bounding each dimension

• vmax (Union[Number, np.ndarray]) – minimum/maximum value bounding each dimension

• perturb_std (Union[Number, np.ndarray]) – standard deviation of the Gaussians used to randomly perturb the Voronoi sites in each dimension.

• n_dimensions (Number, optional) – number of dimensions. None (default) results in a trans-dimensional discretization, with the dimensionality of the parameter space allowed to vary in the range n_dimensions_min-n_dimensions_max

• n_dimensions_min (Number, optional) – minimum and maximum number of dimensions, by default 1 and 10. These parameters are ignored if n_dimensions is not None, i.e. if the discretization is not trans-dimensional

• n_dimensions_max (Number, optional) – minimum and maximum number of dimensions, by default 1 and 10. These parameters are ignored if n_dimensions is not None, i.e. if the discretization is not trans-dimensional

• n_dimensions_init_range (Number, optional) –

  percentage of the range n_dimensions_min` - n_dimensions_max used to initialize the number of dimensions (0.3. by default). For example, if n_dimensions_min = 1, n_dimensions_max = 10, and n_dimensions_init_range = 0.5, the maximum number of dimensions at the initialization is

  int((n_dimensions_max - n_dimensions_min) * n_dimensions_init_range + n_dimensions_max)

• parameters (List[Prior], optional) – a list of free parameters, by default None

• birth_from ({"prior", "neighbour"}, optional) – whether to initialize the free parameters associated with the newborn Voronoi cell by randomly drawing from their prior or by perturbing the value found in the nearest Voronoi cell (default).

Reference Details

is_leaf

name

parameters

all the free parameters defined in this parameter space

perturbation_funcs

the list of perturbation functions allowed in the parameter space linked to the Voronoi discretization. Each function takes in a state (see State) and returns a new state along with the corresponding partial acceptance probability,

\[\begin{split}\underbrace{\alpha_{p}}_{\begin{array}{c} \text{Partial} \\ \text{acceptance} \\ \text{probability} \end{array}} = \underbrace{\frac{p\left({\bf m'}\right)}{p\left({\bf m}\right)}}_{\text{Prior ratio}} \underbrace{\frac{q\left({\bf m} \mid {\bf m'}\right)}{q\left({\bf m'} \mid {\bf m}\right)}}_{\text{Proposal ratio}} \underbrace{\lvert \mathbf{J} \rvert}_{\begin{array}{c} \text{Jacobian} \\ \text{determinant} \end{array}},\end{split}\]

perturbation_weights

a list of perturbation weights, corresponding to each of the perturbation_funcs() that determines the probability of each of them to be chosen during each step

The weights are not normalized and have the following default values:

• Birth/Death perturbations: 1

• Parameter values perturbation: 3

• Voronoi site perturbation: 1

trans_d

indicates whether the current configuration allows changes in dimensionality

add_hyper_params(hyper_params)

Sets the attributes from the given dict and checks for errors

Parameters:

hyper_params (Dict[str, Union[Number, np.ndarray]]) – dictionary of attributes to be set

birth(old_ps_state)

creates a new Voronoi cell, initializes all free parameters associated with it, and returns the pertubed state along with the log of the corresponding partial acceptance probability,

\[\begin{split}\underbrace{\alpha_{p}}_{\begin{array}{c} \text{Partial} \\ \text{acceptance} \\ \text{probability} \end{array}} = \underbrace{\frac{p\left({\bf m'}\right)}{p\left({\bf m}\right)}}_{\text{Prior ratio}} \underbrace{\frac{q\left({\bf m} \mid {\bf m'}\right)}{q\left({\bf m'} \mid {\bf m}\right)}}_{\text{Proposal ratio}} \underbrace{\lvert \mathbf{J} \rvert}_{\begin{array}{c} \text{Jacobian} \\ \text{determinant} \end{array}}.\end{split}\]

In this case, the prior probability of the model \({\bf m}\) is

\[p({\bf m}) = p({\bf c} \mid k) p(k) \prod_i{p({\bf v}_i \mid {\bf c})} ,\]

where \(k\) denotes the number of Voronoi cells, each entry of the vector \({\bf c}\) corresponds to the position of a Voronoi site, and each \(i\)th free parameter \({\bf v}\) has the same dimensionality as \({\bf c}\).

Following [1], \(p({\bf c} \mid k) = \frac{k! \left(N - k \right)!}{N!}\). If we then assume that \(p(k) = \frac{1}{\Delta k}\), where \(\Delta k = k_{max} - k_{min}\), the prior ratio reads

\[\frac{p({\bf m'})}{p({\bf m})} = \frac{(k+1) \prod_i p(v_i^{k+1})}{(N-k)},\]

where \(p(v_i^{k+1})\) denotes the prior probability of the newly born \(i\)th parameter, which may be dependent on \({\bf c}\). The proposal ratio reads

\[\frac{q({\bf m} \mid {\bf m'})}{q({\bf m'} \mid {\bf m})} = \frac{(N-k)}{(k+1) \prod_i q_{v_i}^{k+1}},\]

where \(q_{v_i}^{k+1}\) denotes the proposal probability for the newly born \(i\)th parameter in the new dimension. It is easy to show that, in the case of a birth from neighbor [1] or a birth from prior [2] (see birth_from), \(\lvert \mathbf{J} \rvert = 1\) and \(\alpha_{p} = \frac{p({\bf m'})}{p({\bf m})} \frac{q({\bf m} \mid {\bf m'})}{q({\bf m'} \mid {\bf m})}\). It follows that

\[\alpha_{p} = \frac{(k+1) \prod_i p(v_i^{k+1})}{(N-k)} \frac{(N-k)}{(k+1) \prod_i q_{v_i}^{k+1}} = \frac{\prod_i p(v_i^{k+1})}{\prod_i{q_{v_i}^{k+1}}}.\]

In the case of a birth from prior, \(q_{v_i}^{k+1} = p(v_i^{k+1})\) and

\[\alpha_{p} = \frac{\prod_i p(v_i^{k+1})}{\prod_i{p(v_i^{k+1})}} = 1.\]

In the case of a birth from neighbor, \(q_{v_i}^{k+1} = \frac{1}{\theta \sqrt{2 \pi}} \exp \lbrace -\frac{\left( v_i^{k+1} - v_i \right)^2}{2\theta^2} \rbrace\), where the newly born value, \(v_i^{k+1}\), is generated by perturbing the original value, \(v_i\), of the \(i\)th parameter. This is achieved through a random deviate from the normal distribution \(\mathcal{N}(v_i, \theta)\), with \(\theta\) denoting the standard deviation of the Gaussian used to carry out the perturbation (see, for example, bayesbay.prior.UniformPrior.perturb_std) . The partial acceptance probability is then computed numerically.

Parameters:

old_ps_state (ParameterSpaceState) – current parameter space state

Returns:

• ParameterSpaceState – new parameter space state

• Number – log of the partial acceptance probability, \(log(\alpha_{p}) = \log(\frac{\prod_i p(v_i^{k+1})}{\prod_i{q_{v_i}^{k+1}}})\)

References

[1] (1,2)

Bodin et al. 2012, Transdimensional inversion of receiver functions and surface wave dispersion

[2]

Hawkins and Sambridge 2015, Geophysical imaging using trans-dimensional trees

death(old_ps_state)

removes a new Voronoi cell and returns the pertubed state along with the log of the corresponding partial acceptance probability,

\[\begin{split}\underbrace{\alpha_{p}}_{\begin{array}{c} \text{Partial} \\ \text{acceptance} \\ \text{probability} \end{array}} = \underbrace{\frac{p\left({\bf m'}\right)}{p\left({\bf m}\right)}}_{\text{Prior ratio}} \underbrace{\frac{q\left({\bf m} \mid {\bf m'}\right)}{q\left({\bf m'} \mid {\bf m}\right)}}_{\text{Proposal ratio}} \underbrace{\lvert \mathbf{J} \rvert}_{\begin{array}{c} \text{Jacobian} \\ \text{determinant} \end{array}}.\end{split}\]

It is straightforward to show that this equals the reciprocal of the partial acceptance probability obtained in the case of a birth perturbation (see birth()), i.e.,

\[\alpha_{p} = \frac{\prod_i{q_{v_i}^{k+1}}}{\prod_i p(v_i^{k+1})}.\]

Parameters:

old_ps_state (ParameterSpaceState) – current parameter space state

Returns:

• ParameterSpaceState – new parameter space state

• Number – log of the partial acceptance probability, \(log(\alpha_{p}) = -\log(\frac{\prod_i p(v_i^{k+1})}{\prod_i{q_{v_i}^{k+1}}})\)

get_hyper_param(hyper_param, position=None)

Retrieves the value corresponding to the specified attribute, which may be a function of position

Parameters:

• hyper_param (str) – the name of the attribute

• position (Union[np.ndarray, Number], optional) – the position (in the discretization domain) associated with the value, None by default

Returns:

value corresponding to the specified attribute

Return type:

Union[Number, np.ndarray]

get_perturb_std(*args)

get the standard deviation of the Gaussian used to perturb the discretization

get_perturb_std_birth(position=None)

get the standard deviation of the Gaussian used to perturb the parameter at birth, which may be dependent on the position in the discretization domain

Parameters:

position (Union[Number, np.ndarray], optional) – the position in the discretization domain at which the standard deviation of the Gaussian used to perturb the parameter at birth is returned. Default is None

Returns:

standard deviation of the Gaussian used to perturb the parameter at birth, possibly at the specified position

Return type:

Number

has_hyper_param(hyper_param)

Whether or not the Prior instance has the specified attribute

Parameters:

hyper_param (str) – the name of the attribute

initialize(position=None)

initializes the parameter space linked to the Voronoi tessellation

Returns:

an initial parameter space state, or a list of parameter space states

Return type:

Union[ParameterSpaceState, List[ParameterSpaceState]

log_prior(*args)

BayesBay implements the grid trick, which calculates the prior probability of a Voronoi discretization through the combinatorial formula \({N \choose k}^{-1}\), with k denoting the number of Voronoi sites and N the number of possible positions allowed for the sites [3].

References

[3]

Bodin and Sambridge (2009), Seismic tomography with the reversible jump algorithm

log_prob_initialize_discretization(ps_state)

The log of the partial acceptance probability of the birth of the discretization. This includes only the discretization but not the parameter values.

\(\frac{p(k')}{p(k)} \frac{p(c'|k')}{p(c|k)} \frac{q(c'|m)}{q(c|m')}\)

Parameters:

ps_state (ParameterSpaceState) – the newly-born parameter space state

Returns:

the log of the partial acceptance probability of the birth of the discretization

Return type:

Number

nearest_neighbour(discretization, query_point)

returns the index of the nearest neighbour of a given query point in the discretization

Parameters:

• discretization (np.ndarray) – the discretization

• query_point (Union[Number, np.ndarray]) – the query point

Returns:

the index of the nearest neighbour point

Return type:

int

perturb_value(old_ps_state, isite)

perturbs the value of one Voronoi site and calculates the log of the partial acceptance probability

\[\begin{split}\underbrace{\alpha_{p}}_{\begin{array}{c} \text{Partial} \\ \text{acceptance} \\ \text{probability} \end{array}} = \underbrace{\frac{p\left({\bf m'}\right)}{p\left({\bf m}\right)}}_{\text{Prior ratio}} \underbrace{\frac{q\left({\bf m} \mid {\bf m'}\right)}{q\left({\bf m'} \mid {\bf m}\right)}}_{\text{Proposal ratio}} \underbrace{\lvert \mathbf{J} \rvert}_{\begin{array}{c} \text{Jacobian} \\ \text{determinant} \end{array}}.\end{split}\]

Parameters:

• old_ps_state (ParameterSpaceState) – the current parameter space state

• isite (Number) – the index of the Voronoi site to be perturbed

Returns:

the new parameter space state and its associated partial acceptance probability excluding log likelihood ratio

Return type:

Tuple[ParameterSpaceState, Number]

sample(*args)

sample a random ParameterSpaceState instance, including the number of dimensions, the Voronoi sites, and the parameter values

sample_discretization()

sample a parameter space state that only contains the discretization

sample_site()

draws a Voronoi-site position at random within the discretization domain

set_custom_initialize(initialize_func)

sets a custom initialization function

Parameters:

initialize_func (Callable[["Prior", np.ndarray], np.ndarray]) – The function to use for initialization. This function should take a Prior instance and optionally an array of positions as input arguments, and produce an array of values as output.

Examples

def my_init(
    param: bb.prior.Prior,
    position: np.ndarray
) -> np.ndarray:
    print("This is my custom init!")
    return np.ones(len(position))

my_param.set_custom_initialize(my_init)

back to top