API reference

Bayesian MSD fitting

Any valid MSD function MSD(Δt) defines a stationary Gaussian process with the appropriate correlation structure. Therefore, instead of fitting the graph of “empirically” calculated MSDs, we can perform Bayesian inference of parametric MSD curves, using the likelihood function defined through these Gaussian processes.

We discriminate two cases, depending on what exactly is stationary:

  • the trajectories themselves might be sampled from a stationary process (e.g. distance of two points on a polymer). In terms of the MSD, the decisive criterion is MSD(inf) < inf. In this case, it is straightforward to prove the following relation between the MSD μ and autocovariance γ of the process:

    μ(k) = 2*( γ(0) - γ(k) )

    Thus, the full autocovariance function can be obtained from the MSD and the steady state covariance γ(0). For decaying correlations (γ(k) --> 0 as k --> ∞) we furthermore see that 2*γ(0) = μ(∞) is the asymptotic value of the MSD. Finally, this allows us to calculate the covariance matrix of the process as

    C_ij := <x_i*x_j>
      = γ(|i-j|)
      = γ(0) - 1/2*μ(|i-j|)

    We call this case a steady state of order 0.

  • a special case of order 0 is the scenario where we have a stationary process, but when evaluating e.g. likelihoods we wish to condition on the first data point of the trajectory. This renders the resulting likelihoods comparable to those of an order 1 stationary process and can thus be useful sometimes; we refer to this scenario as ss_order = 0.5.

  • in many cases (e.g. sampling a diffusing particle’s position) the trajectories themselves will not be stationary, but the increment process is. In this case the Gaussian process of interest is the one generating the increments of the trajectory, whose autocorrelation is the second derivative of the MSD:

    γ(k) = 1/2 * (d/dk)^2 μ(k)

    where derivatives should be understood in a weak (“distributional”) sense. More straightforwardly, the correlation matrix of the increments is given by

    C_ij := <(x_{i+1} - x_i)(x_{j+1}-x_j)>
      = 1/2 * ( μ(t_{i+1} - t_j) + μ(t_i - t_{j+1})
               -μ(t_{i+1} - t_{j+1}) - μ(t_i - t_j) )

    where by definition we let μ(-k) = μ(k). In this case, we talk about a steady state of order 1.

In any case, the covariance matrix C (potentially together with a mean / drift term for steady state order 0 / 1 respectively) defines a Gaussian process, which lets us assign a likelihood to the generating MSD. Via this construction, we can perform rigorous Bayesian analysis, cast in the familiar language of MSDs.

This package provides a base class for performing such inferences / fits, namely Fit. We also provide a few example implementations of specific fitting schemes in the lib submodule. Finally, the Profiler allows to explore the posterior once a point estimate has been found, by tracing out either conditional posterior or profile posterior curves in each parameter direction. Note that you can also just sample the posterior by MCMC.

See also

Fit, Profiler

bayesmsd.deco

Decorators to use when defining custom MSD functions

There is some functionality that should be common to all/most functions to be used in Fit.params2msdm implementations. This module provides decorators implementing this common functionality. Specifically, these fall into two categories:

  • technical convenience, i.e. making sure that the argument is properly cast to a numpy array and dealing with dt <= 0 (i.e. MSD(0) == 0 should be true for any MSD function for technical reasons, even if lim_{dt-->0} MSD(dt) != 0). This is taken care of by the MSDfun decorator.

  • contributions to realistic MSD curves due to common imaging artifacts, i.e. localization error and motion blur. These are implemented in the imaging decorator.

Taking both together, the definition for a simple powerlaw MSD function would take the following form:

>>> G, alpha = ... # from input
...
... @bayesmsd.deco.MSDfun
... @bayesmsd.deco.imaging(noise2=..., f=..., alpha0=alpha)
... def powerlawMSD(dt, G=G, alpha=alpha):
...     return G*dt**alpha

Note that MSDfun should always be the topmost decorator, since it handles user input. The imaging decorator first adds motion blur—parametrized by the fractional exposure f which is given experimentally and the early time (true) MSD scaling alpha0, which should be determined from the true MSD curve—and then adds the (squared) localization error noise2. Overall, there is thus just one additional free parameter, namely the localization error.

Note how we pass G and alpha (which are defined in the surrounding namespace) to the function as default arguments. This gets the scoping correct, since default arguments are evaluated at definition time (as opposed to just using the externally defined variables, whose values might change between definition and execution of the function).

See also

bayesmsd

bayesmsd.deco.method_verbosity_patch(meth)

(internal) Decorator for class methods, temporarily changing verbosity

bayesmsd.deco.MSDfun(fun)

Decorator for MSD functions

This is a decorator to use when implementing Fit.params2msdm. It takes over some of the generic polishing: it assumes that the decorated function has the signature function(np.array) --> np.array and

  • ensures that the argument is cast to an array if necessary (such that you can then also call msd(5) instead of msd(np.array([5]))

  • ensures that dt > 0 by taking an absolute value and setting msd[dt==0] = 0 without calling the wrapped function. You can thus ignore the dt == 0 case in implementing an MSD function. Note that msd[dt==0] = 0 should always be true. However, e.g. in the case of localization error, we might have lim_{Δt-->0} MSD(Δt) = 2σ² != 0.

See also

Fit, imaging

bayesmsd.deco.imaging(noise2=0, f=0, alpha0=1)

Add imaging artifacts (localization error & motion blur) to MSDs.

This decorator should be used when defining MSD functions for SPT experiments, to add artifacts due to the imaging process. These are a) localization error and b) motion blur, caused by finite exposure times. Use this decorator after MSDfun, like so:

>>> @bayesmsd.deco.MSDfun
... @bayesmsd.deco.imaging(...)
... def msd(dt):
...     ...
Parameters:

  • noise2 (float >= 0) – the variance (σ²) of the Gaussian localization error to add

  • f (float, 0 <= f <= 1) – the exposure time as fraction of the frame time.

  • alpha0 (float, 0 <= alpha0 <= 2 or 'auto') – the effective short time scaling exponent. Set to 'auto' to determine from finite differences of the “raw” MSD around Δt=f.

Notes

f = 0 is ideal stroboscopic illumination, i.e. no motion blur. Accordingly, the value of alpha0 is not used in this case.

See also

MSDfun

bayesmsd.fit

Implementation of Fit base class

See also

bayesmsd, Fit, Profiler

class bayesmsd.fit.Fit(data)

Bases: object

Abstract base class for MSD fits; the backbone of bayesmsd.

Subclass this to implement fitting of a specific functional form of the MSD. See also the existing library of fits in bayesmsd.lib.

Parameters:

data (noctiluca.TaggedSet, pandas.DataFrame, list of numpy.ndarray) –

the data to run the fit on.

This input is handled by the userinput.make_TaggedSet() function of the noctiluca package and thus accepts a range of formats. A dataset of trajectories with N loci, T frames, and d spatial dimensions can be specified as list of numpy arrays with dimensions (N, T, d), (T, d), (T,), where T can vary between trajectories but N and d should be the same. A pandas dataframe should have columns (particle, frame, x1, y1, z1, x2, y2, z2, ...), where particle identifies which trajectory each entry belongs to, which frame is the frame number in which it was detected. For precise specs see noctiluca.util.userinput.

data

the data to use. Note that the current selection in the data is saved internally

Type:

TaggedSet of Trajectory

d

spatial dimension of trajectories in the data

Type:

int

T

maximum length of the trajectories, in frames

Type:

int

ss_order

steady state order. Often this will be a fixed constant for a particular Fit subclass.

Type:

{0, 1}

parameters

the parameters for this fit. Each entry should carry a sensible name and be an instance of bayesmsd.Parameter. By default, this is populated with a parameter m1 for each dimension, which is fixed to 0. This can be used for the first moment (mean/drift) in the fits, or removed/overridden.

Type:

dict

constraints

allows to specify constraints on the parameters that will be implemented as smooth penalty on the likelihood. Can also take care of feasibility constraints: by default, there is a constraint checking that the covariance matrix given by the current MSD is positive definite (otherwise we could not even evaluate the likelihood). See Notes section.

Type:

list of constraint functions

max_penalty

constant cutoff for the penalty mechanism. Infeasible sets of parameters (where the likelihood function is not even well-defined) will be penalized with this value. For any set of parameters penalized with at least this value, likelihood evaluation is skipped and the value max_penalty assigned as value of the minimization target. Default value is 1e10, there should be little reason to change this.

Type:

float

verbosity

controls output during fitting. 0: no output; 1: error messages only; 2: informational; 3: debugging

Type:

{0, 1, 2, 3}

Notes

constraints are functions with signature constraint(params) --> float. The output is interpreted as follows:

  • x <= 0 : infeasible; maximum penalization

  • 0 <= x <= 1 : smooth penalization: penalty = exp(cot(pi*x))

  • 1 <= x : feasible; no penalization

Thus, if e.g. some some function fun of the parameters should be constrained to be positive, you would use fun(params)/eps as the constraint, with eps some small value setting the tolerance region. If there are multiple constraints, always the strongest one is used. For infeasible parameters, the likelihood function is not evaluated; instead the “likelihood” is just set to -Fit.max_penalty.

Note that there is a default constraint checking positivity of the covariance matrix. If your functional form of the MSD is guaranteed to satisfy this (e.g. for a physical model), you can remove this constraint for performance.

Upon subclassing, it is highly recommended to initialize the base class first thing:

>>> def SomeFit(Fit):
...     def __init__(self, data, *other_args):
...         super().__init__(data) # <--- don't forget!

This class uses scipy.optimize.minimize to find the MAP parameter estimate (or MLE if you leave logprior flat). When running the fit, you can choose between the simplex (Nelder-Mead) algorithm or gradient descent (L-BFGS-B). The latter uses the stopping criterion f^k - f^{k+1}/max{|f^k|,|f^{k+1}|,1} <= ftol, which is inappropriate for log-likelihoods (which should be optimized to fixed accuracy of O(0.1) independent of the absolute value, which in turn might be very large). We therefore use Nelder-Mead by default, which does not depend on derivatives and thus also has an absolute stopping criterion.

If this function runs close to a maximum, e.g. in Profiler or when using successive optimization steps, we can fix the problem with gradient-based optimization by removing the large absolute value offset from the log-likelihood. This functionality is also exposed to the end user, who can overwrite the initial_offset method to give a non-zero offset together with the initial values provided via initial_params.

The rationale for pre-populating parameters with the first moment parameters is that these are rarely interesting when writing a new fit, so wouldn’t be implemented by a lazy developer (yours truly) by default. Just carrying them through, however, is completely trivial and exposes the trend fitting functionality to the end user by simply “unfixing” the corresponding parameters (i.e. setting .fix_to = None).

vprint(v, *args, **kwargs)

Prints only if self.verbosity >= v.

abstract params2msdm(params)

Definition of MSD in terms of parameters

This is the core of the fit definition. It should give a list of tuples as required by GP.ds_logL

Parameters:

params (dict) – the current parameter values

Returns:

list of tuples (msd, m) – for the definition of msd, use of the MSDfun decorator (and potentially also imaging) is recommended.

See also

GP.ds_logL, MSDfun, imaging

logprior(params)

Prior over the parameters

Use this function if you want to specify a prior over the parameters that will be added to the log-likelihood. Default is a flat prior, i.e. return 0.

Parameters:

params (dict) – the current parameter values; note that only the free parameters are given

Returns:

float

abstract initial_params()

Give initial values for the parameters

You can use self.data to perform some ad hoc estimation (e.g. from the empirical MSD, using noctiluca.analysis.MSD) or just return constants.

Returns:

params (dict) – initial values for the parameters. Should have the same keys as self.parameters.

initial_offset()

Log-likelihood offset associated with initial parameters

See Notes section of class documentation.

Returns:

float

See also

Fit

constraint_Cpositive(params)

Constraint for positive definiteness of covariance matrix

This can serve as an example of a non-trivial constraint. Note that you may not want to use it if positivity is already guaranteed by the functional form of your MSD. See also Notes section of class doc.

Returns:

float

Notes

This function checks whether the spatial components are identical using python’s is. So if you are implementing an MSD with identical spatial components, you should return the final list as self.d*[(msd, mean)], such that this constraint checks positivity only once.

See also

Fit

MSD(params, dt=None)

Return the (trend-free) MSD evaluated at dt. Convenience function.

Parameters:

  • params (dict) – the current parameter values. Like the 'params' entry of the output dict from any fit.

  • dt (array-like) – the time lags (in frames) at which to evaluate the MSD. If left unspecified, we return a callable MSD function

  • kwargs (additional keyword arguments) – are forwarded to the MSD functions returned by params2msdm.

Returns:

callable or np.array – the MSD function (summed over all dimensions), evaluated at dt if provided.

Notes

This function returns the “trend-free” MSD; meaning we ignore potential values for the mean/trend parameters m1. For ss_order == 0 this is correct behavior either way; for ss_order == 1 the trend term can be added as (m1*Δt)^2.

expand_fix_values(fix_values=None)

Assemble full list of values to fix (internal + given)

Merge the given fix_values with the parameter-level specifications and ensure that all entries are appropriate callables (can also be specified as numerical constants or other parameter names)

Parameters:

fix_values (dict, optional) – values to fix, beyond what’s already in self.parameters[...].fix_to

Returns:

fix_values (dict) – same as input, plus internals, and constants resolved

Notes

Parameter resolution order:

  • (marginalization)

  • free parameters

  • fixed to constant

  • fixed to other parameter by name

  • fixed to callable

So you cannot fix to another parameter that is itself fixed to a callable. You can use parameter.fix_to = other_parameter.fix_to instead to just copy the function over instead of fixing by name. This helps prevent undetectable cycles in the fixes.

independent_parameters(fix_values=None)

Return the names of the independent (not fixed) parameters

Parameters:

fix_values (dict) – should be the same as was / will be handed to run

Returns:

list of str – a list with the names of the independent parameters, i.e. those that are not fixed to some other value

class MinTarget(fit, *, fix_values=None, offset=0, adjust_prior_for_fixed_values=True)

Bases: object

Helper class; acts as cost function for optimization.

Beyond acting as the actual cost function through implementation of the () operator, this class also defines the conversion between parameter dicts (where parameters are named and unordered) and parameter arrays (which are used for optimization). Most of this is for internal use.

Parameters:

  • fit (Fit) – the Fit object that this target is associated with.

  • fix_values (dict) – any additional fixes for parameters, beyond what’s already in fit.parameters[...].fix_to.

  • offset (float) – the global offset to subtract; see Notes section of Fit.

  • adjust_prior_for_fixed_values (bool) – if False, pretend that the parameters in fix_values are still free. Important for the Profiler.

fit

the Fit object that this target is associated with.

Type:

Fit

fix_values

the complete set of values to keep fixed; keys are parameter names, values are (guaranteed to be) callables.

Type:

dict

fix_values_resolution_order

list of parameter names. The fixes in fix_values will be executed in this order.

Type:

list of str

param_names

list of parameter names that remain independent; this determines the order of entries when converting parameter dicts to arrays.

Type:

list of str

offset

constant to be subtracted from the target (i.e. added to the likelihood). Can be used to ensure that the minimum value of the target is close to 0.

Type:

float

params_array2dict(params_array)

Convert a parameter array to dict

Parameters:

params_array (np.ndarray)

Returns:

dict

See also

params_dict2array

params_dict2array(params)

Convert a parameter dict to array

Parameters:

params_dict (dict)

Returns:

np.ndarray

See also

params_array2dict

__call__(params_array)

Actual minimzation target

This function evaluates penalty, prior, and likelihood for the current parameters and returns the negative of their sum. If penalty < 0, instead of evaluating the likelihood, we just return the maximum penalization.

Parameters:

params_array (np.ndarray)

Returns:

float

profile_marginalized_params(params)

Run the profiler over the marginalized parameters

Parameters:

params (dict)

Returns:

dict – for each marginalized parameter: point estimate and 95% credible interval

run(init_from=None, optimization_steps=('simplex',), maxfev=None, fix_values=None, adjust_prior_for_fixed_values=True, full_output=False, show_progress=False)

Run the fit

Parameters:

  • init_from (dict) – initial point for the fit, as a dict with fields 'params' and 'logL', like the ones this function returns. The 'logL' entry is optional; you can also leave it out or set to 0.

  • optimization_steps (tuple) – successive optimization steps to perform. Entries should be 'simplex' for Nelder-Mead, 'gradient' for gradient descent, or a dict whose entries will be passed to scipy.optimize.minimize as keyword arguments.

  • maxfev (int or None) – limit on function evaluations for 'simplex' or 'gradient' optimization steps

  • fix_values (dict) – can be used to keep some parameter values fixed or express them as function of the other parameters. See class doc for more details.

  • adjust_prior_for_fixed_values (bool) – whether to evaluate the model prior for the restricted parameter space resulting from fix_values or over the full parameter space. While the former is “correct” in most cases, the latter is important if we want to explore the likelihood landscape around a previously found optimum, as the Profiler does.

  • full_output (bool) – Set to True to return the output dict (c.f. Returns) and the full output from scipy.optimize.minimize for each optimization step. Otherwise (full_output == False, the default) only the output dict from the final run is returned.

  • show_progress (bool) – display a tqdm progress bar while fitting

  • verbosity ({None, 0, 1, 2, 3}) – if not None, overwrites the internal self.verbosity for this run. Use to silence or get more details of what’s happening

Returns:

dict – with fields 'params', a complete set of parameters; 'logL', the associated value of the likelihood (or posterior, if the prior is non-trivial).

evidence(show_progress=False, conf=0.8, conf_tol=0.1, n_cred=10, n_steps_per_cred=2, f_integrate=0.99, log10L_improper=3)

Estimate evidence for this model

The parameters to this function are mostly technical and can remain at their default values for most use cases. It is possible that in some cases speedups can be achieved by tuning them, but generally the defaults should be reasonable enough.

Parameters:

  • show_progress (bool) – display progress bar(s)

  • conf (float) – confidence level for the initial Profiler run. Will usually be some not-too-high value.

  • conf_tol (float) – tolerance for conf. Usually relatively high.

  • n_cred (float) – how far (at most) from the point estimate to evaluate the likelihood function, in multiples of the initial Profiler credible interval.

  • n_steps_per_cred (float) – how many grid points to place between point estimate and boundary of initial credible interval.

  • f_integrate (float) – proxy parameter to determine the threshold until which the likelihood will be evaluated. Can be interpreted as “if the posterior was Gaussian, we would cover at least a fraction f of its mass”.

  • log10L_improper (float) – width of surrogate prior for improper priors. See Notes.

Returns:

float – estimated evidence for this model

Notes

This function estimates model evidence by integrating the likelihood over a non-uniform grid. The outline of the algorithm is as follows:

  • run the fit to get a point estimate

  • run a Profiler (in conditional posterior mode) to get order of magnitude estimates for the steps in different directions over which the likelihood changes. This run is controlled by conf and conf_tol; since we do not need precise estimates, but just an idea of how big steps should generally be, the default is to search with relatively low accuracy (high conf_tol) for a not-too-strong decrease in likelihood (medium values for conf).

  • given the profiler results, assemble a grid of parameter values over which the likelihood may be evaluated. This step is controlled by n_cred (how big to make the grid) and n_steps_per_cred (how fine to make the grid)

  • start by evaluating the grid completely within the bounds given by the profiler

  • grow the integration region iteratively: at each step, find the grid points that are adjacent to another one that a) has been evaluated already and b) has a value above the cutoff controlled by f_integrate. Evaluate the likelihood on all these candidate points; repeat

  • with the likelihood evaluated at all relevant grid points, integrate numerically by multiplying the voxel volume and summing.

For the resulting evidence value to be accurate and comparable across different models, it is tantamount that the priors be properly normalized. This means that a) use this function only with models that implement the logprior function correctly; b) models with improper priors (e.g. for localization error we usually use a log-flat prior) are problematic. We sneak around the latter issue with a debatable method: for the cases where we have improper priors, we fix a suitable, broad, and importantly proper prior on the location of the point estimate. This approach is clearly suspicious, because we fix the prior to a result from the data. One might argue that at the end of the day we are basically “just” fixing the order of magnitude (think: units) of a given parameter, which we might also do by checking the experimental protocol. But at the end of the day, it remains a questionable method; any reader with better ideas for how to handle this case should email me.

Having accepted the dubious approach to improper priors: the surrogate (proper) prior we finally use is a Gaussian with standard deviation log(10)*log10L_improper. The idea here is that if the improper prior is log-flat (my most common use case), the parameter log10L_improper basically gives the number of valid digits in a parameter estimate that we would accept as “not fine tuned”. Specifically: consider “fitting” a model without any parameters; if I can improve upon that fit by introducing one new parameter with the value 17, we might accept the second model as reasonable; if the fit only improves if I fix the new parameter to 17.2193427, we might say that the new model requires too much fine tuning. Upon personal reflection, log10L_improper = 3 valid digits seems to be reasonable; reasonable people might differ.

bayesmsd.fitgroup

class bayesmsd.fitgroup.FitGroup(fits_dict)

Bases: Fit

Running interdependent fits

The FitGroup allows running fits to different data sets simultaneously, while fixing parameters across these fits. A common use case is a scenario where we acquired SPT data of the same biological system with different technical settings: choice of microscope, choice of lag time, localization error & motion blur, etc. In this case we would like to run a fit that fits the same “true” MSD to the data, while maintaining independence for e.g. the localization errors associated with each data set.

Usage is (hopefully) relatively intuitive:

  • define a few “simple” fits and collect them in a dict

  • give that dict to FitGroup

  • use the usual .fix_to = ... mechanism to fix parameters in the fit group. All parameters from the individual fits can be used, with their names prefixed by the fit name (the key for that fit in the dict)

  • FitGroup.run() as you would Fit.run()

Parameters:

fits_dict (dict of Fit) – the fits to run; we refer to the dict keys as “fit names”

parameters

collection of all the individual fit parameters, with their names prefixed by the fit name (so 'α (dim 0)' of the fit 'fit1' would become 'fit1 α (dim 0)')

Type:

dict of Parameter

constraints

same as Fit.constraints. Note that you can use constraints within the individual fits, as well as on the group level. It usually makes sense to stay “as local as possible”, though this is not a requirement.

Type:

list of callables

max_penalty

same as Fit.max_penalty

Type:

float

verbosity

same as Fit.verbosity

Type:

int

improper_priors

aggregated list of Fit.improper_priors, with names adjusted to the joint fit (i.e. prefixed with the fit name)

Type:

list of str

joint_param_prefix(fitname)
make_joint_param_name(fitname, paramname)
make_fit_param_name(fitname, paramname)
params2msdm(params)

Definition of MSD in terms of parameters

This is the core of the fit definition. It should give a list of tuples as required by GP.ds_logL

Parameters:

params (dict) – the current parameter values

Returns:

list of tuples (msd, m) – for the definition of msd, use of the MSDfun decorator (and potentially also imaging) is recommended.

See also

GP.ds_logL, MSDfun, imaging

MSD(params, dt)

Return the (trend-free) MSD evaluated at dt. Convenience function.

Parameters:

  • params (dict) – the current parameter values. Like the 'params' entry of the output dict from any fit.

  • dt (array-like) – the time lags (in frames) at which to evaluate the MSD. If left unspecified, we return a callable MSD function

  • kwargs (additional keyword arguments) – are forwarded to the MSD functions returned by params2msdm.

Returns:

callable or np.array – the MSD function (summed over all dimensions), evaluated at dt if provided.

Notes

This function returns the “trend-free” MSD; meaning we ignore potential values for the mean/trend parameters m1. For ss_order == 0 this is correct behavior either way; for ss_order == 1 the trend term can be added as (m1*Δt)^2.

initial_params()

Aggregate initial parameters from individual fits

Returns:

dict – initial parameters from individual fits, with properly prefixed names

logprior(params)

Aggregate log-priors from individual fits

Parameters:

params (dict) – parameters for which to evaluate the prior

Returns:

float – aggregated log-prior

Notes

Which parameters are handed in with params determines what is considered a free parameter; so ensure to not hand in parameters that factually are fixed.

The ability of a FitGroup to calculate its prior is not really necessary and not used in the code; this function is thus implemented mostly for completeness/rigor (it is not necessary for Fit.evidence() to work; what we need there is just that the individual fits take care of their respective priors respectively).

class MinTarget(*args, **kwargs)

Bases: MinTarget

Cost function for the FitGroup; overrides Fit.MinTarget

Mainly take care of proper parameter and prior management; actual likelihood evaluations are delegated to the individual fits.

Parameters are same as for Fit.MinTarget.

bayesmsd.gp

Implementation of Gaussian process logic

This module is mostly for internal use of the bayesmsd package. It covers the definition of (stationary) Gaussian processes in terms of their MSD, and the associated likelihood function ds_logL.

Gaussian processes can also serve as generative model, i.e. one can sample trajectories with a given MSD curve. This is implemented in the generate function in this module.

See also

bayesmsd, ds_logL, generate, generate_dataset_like

bayesmsd.gp.msd2C_fun(msd, ti, ss_order)

msd2C for MSDs expressed as functions / non-integer times.

Parameters:

  • msd (callable, use the MSDfun decorator) – note that for ss_order < 1 we expect msd(np.inf) to be well-defined.

  • ti (np.ndarray, dtype=int) – times at which there are data in the trajectory

  • ss_order ({0, 1}) – steady state order. See module documentation.

Returns:

np.ndarray – covariance matrix

See also

MSDfun, bayesmsd

class bayesmsd.gp.GP

Bases: object

exception BadCovarianceError

Bases: RuntimeError

logL(ss_order, msd, mean=0)

Gaussian process likelihood for a given trace

Parameters:

  • trace ((T,) np.ndarray) – the data. Should be recorded at constant time lag; missing data are indicated with np.nan.

  • ss_order ({0, 0.5, 1}) – steady state order; see module documentation.

  • msd (np.ndarray) – the MSD defining the Gaussian process, evaluated up to (at least) T = len(trace).

  • mean (float) – the first moment of the Gaussian process. For ss_order == 0 this is the mean (i.e. same units as the trajectory), for ss_order == 1 this is the mean of the increment process, i.e. has units trajectory per time.

Returns:

float – the log-likelihood

See also

ds_logL

ds_logL(ss_order, msd_ms)

Gaussian process likelihood on a data set.

Parameters:

  • data (a TaggedSet of Trajectory)

  • ss_order ({0, 0.5, 1}) – steady state order; see module documentation.

  • msd_ms (list of tuples (msd, mean)) – this should be a list with one entry for each spatial dimension of the data set. The first entry of each tuple is a function, ideally decorated with MSDfun. The second should be a float giving the mean/drift for the process (often zero).

Returns:

float – the total log-likelihood of the data under the Gaussian process specified by the given MSD.

Notes

Parallel-aware (unordered). However, in practice I find that it is usually faster to parallelize runs over multiple data sets / with different parameters, if possible. In a benchmarking run for the internal parallelization in this function I saw no observable benefit of running on more than 5 cores (presumably due to overhead in moving data to the workers).

See also

Fit, noctiluca.parallel

bayesmsd.gp.generate(msd_def, T, n=1)

Sample trajectories from a given MSD / fitparameters

Parameters:

  • msd_def (tuple) –

    one of two options:

    • (fit, res) where fit is an instance of Fit (e.g. one of the classes from lib) and res is the result of a previous run of fit (i.e. a dict with 'params' and 'logL' entries; the latter is not used here)

    • (msdfun, ss_order, d) where msdfun is a callable (ideally wrapped with the MSDfun decorator), ss_order in {0, 1} (see module doc), and d is the spatial dimension of the trajectories

  • T (int) – the length of the trajectories to be generated, in frames

  • n (int, optional) – the number of trajectories to generate

Returns:

TaggedSet – a data set containing the drawn trajectories

bayesmsd.gp.generate_dataset_like(data, msd_def)

Create a dataset like the reference, but from given MSD.

The returned TaggedSet contains trajectories that match the input dataset in number, length, and missing frames, but are sampled from the MSD defined by msd_def (which is usually a fit to these data).

Parameters:

  • data (TaggedSet of Trajectories) – the input dataset. Note that the actual data will not be used, just meta info like number and length of trajectories, and missing frames

  • msd_def ((fit, res) or (msd_fun, ss_order, d)) – the definition of the MSD to sample from; see generate.

Returns:

TaggedSet of Trajectories

Notes

This function assumes N = 1 (single particle) for all trajectories; d (number of spatial dimensions) should be consistent across the dataset.

bayesmsd.lib

A library of Fit implementations

This module provides some commonly used MSD fits. Most notably:

  • Fitting a powerlaw through NPXFit: NPX stands for Noise + Powerlaw + X, i.e. this fit class has the capability of fitting a powerlaw MSD to your data and supplement it with standard imaging artifacts (localization error and motion blur) as well as a freeform extension at long times (we use cubic splines). The latter is disabled (n = 0) by default, thus providing a standard powerlaw fit.

  • Fitting cubic splines to your MSD curve via SplineFit. Doing this for varying number n of spline points allows to discern which features of an empirical MSD curve are significant and which are just noise. To compare fits with different number of spline points, you can use AIC.

class bayesmsd.lib.SplineFit(data, ss_order, n, previous_spline_fit_and_result=None, motion_blur_f=0)

Bases: Fit

Fit a spline MSD

The MSD in this case is parametrized by the positions of a few spline points, between which we interpolate with cubic splines. The boundary conditions for the splines are set such that the fitted MSD extrapolates beyond the data as a powerlaw, except for large times in the ss_order == 0 case, in which the MSD converges to a constant at infinite time. To achieve this within the spline fits, the time coordinate is compactified from t in [1, T] to x in [0, 1]:

  • if ss_order < 1, we need t = inf to be accessible. We therefore choose the compactification x = 4/π*arctan(log(t)/log(T)), such that x(t=∞) = 2 while x(t=T) = 1.

  • if ss_order == 1, we simply work in log-space, and normalize: x = log(t)/log(T). Thus x(t=T) = 1.

For the y-coordinate of the spline we apply a simple log-transform, y = log(MSD). We then use the boundary condition that the second derivative of the spline vanishes, meaning the fitted MSDs can be extrapolated quite naturally by powerlaws.

Parameters:

  • data (noctiluca.TaggedSet, pandas.DataFrame, list of numpy.ndarray) – the data to fit. See Fit.

  • ss_order ({0, 1}) – the steady state order to assume

  • n (int >= 2) – the number of spline points

  • previous_spline_fit_and_result (tuple (SplineFit, dict)) – this can be used to initialize the current fit from the resulting spline of a previous one. Very useful when running the above-mentioned model selection task over spline knots. First entry should be the Fit object used, second one should be the resulting dict with keys 'params' and 'logL'.

  • motion_blur_f (float in [0, 1]) – fractional exposure time (i.e. exposure time as fraction of the time between frames). The resulting motion blur will be taken into account in the fit.

motion_blur_f
Type:

float

n
Type:

int

ss_order
Type:

{0, 1}

logT

log of trajectory length; used in compactification of spline coordinates

Type:

float

upper_bc_type

the boundary condition for the upper end of the spline. For ss_order == 0, this sets the derivative to zero (MSD should plateau at infinity), for ss_order == 1 we use natural boundary conditions, i.e. vanishing second derivative.

Type:

str or tuple

x_last

the theoretical bound on the x coordinate for the spline points; 2 for ss_order < 1, 1 for ss_order == 1.

Type:

{1, 2}

x_max

the maximum value for x such that the corresponding dt is still numerically representable

Type:

float

prev_fit

same as parameter previous_spline_fit_and_result.

Type:

tuple

Notes

Clearly the number of spline points controls the goodness of fit and the degree of overfitting. It is thus recommended to use some information criterion (e.g. AIC) to determine a reasonable level of detail given the data. This approach turns out to be pretty useful in understanding which features of an “empirical MSD” (like calculated by noctiluca.analysis.MSD) are reliable, and which ones are just noise.

The parameters for this fit are the coordinates of the spline knots, in the transformed coordinate system ((dt, MSD) --> (x, y) == (compactify(dt), log(MSD))). Since the total extent of the data along the time axis is fixed, we also fix the first and last x-coordinates, such that ultimately the free parameters are x1 through x(n-2) and y0 through y(n-1)

See also

Fit, bayesmsd

logprior(params)

Prior over the parameters

Use this function if you want to specify a prior over the parameters that will be added to the log-likelihood. Default is a flat prior, i.e. return 0.

Parameters:

params (dict) – the current parameter values; note that only the free parameters are given

Returns:

float

compactify(dt)

Compactification used for the current fit

Parameters:

dt (np.array, dtype=float) – the time lags to calculate compactification for. Should be > 0 but might include np.inf.

See also

SplineFit, decompactify_log

decompactify_log(x)

Decompactify spline points (convenience function)

Parameters:

x (np.array, dtype=float) – the compactified x-coordinates

Returns:

log_dt (np.array) – the corresponding log(dt [frames]) values

See also

compactify

params2msdm(params)

Calculate the current spline MSD

See also

Fit.params2msdm

initial_params()

Give suitable initial parameters for the spline

To find proper initial parameters, we perform a simple powerlaw fit to the empirical MSD. In the ss_order < 1 case this is just used as boundary condition for a two-point spline between the first time lag and infinity (where we use the empirical steady state variance as initial value). If ss_order == 1 this fitted powerlaw is the initial MSD.

Returns:

params (np.ndarray, dtype=float) – the inital spline knots, in the internal x-y-coordinates (i.e. compactified)

See also

Fit.initial_params

initial_offset()

Used when starting from a previous SplineFit

See also

Fit.initial_offset

constraint_dx(params)

Make sure the spline points are properly ordered in x

We impose this constraint mainly to avoid crossing of spline points, which usually leads to the spline diverging. On top of that, conceptually this makes the solution well-defined.

Parameters:

params (np.ndarray, dtype=float) – the current fit parameters

Returns:

float – the constraint score

See also

Fit

constraint_logmsd(params)

Make sure the Spline does not diverge

Parameters:

params (np.ndarray, dtype=float) – the current fit parameters

Returns:

float – the constraint score

See also

Fit

class bayesmsd.lib.NPXFit(data, ss_order, n=0, previous_NPXFit_and_result=None, motion_blur_f=0)

Bases: Fit

(N)oise + (P)owerlaw + (X) (i.e. spline) fit

This scheme can be used to fit powerlaw MSDs with a few extensions: we include localization error and motion blur, and towards long time the MSD might deviate from a powerlaw, which we then fit with a cubic spline (see also SplineFit).

Parameters:

  • data (noctiluca.TaggedSet, pandas.DataFrame, list of numpy.ndarray) – the data to fit. See Fit.

  • ss_order ({0, 1}) – the steady state order to assume

  • n (int >= 0) – the number of spline points; set n = 0 to fit a pure powerlaw.

  • previous_NPXFit_and_result (tuple (NPXFit, dict)) – see also SplineFit; this can be used to run model selection over n.

  • motion_blur_f (float in [0, 1]) – fractional exposure time (i.e. exposure time as fraction of the time between frames). The resulting motion blur will be taken into account in the fit.

motion_blur_f
Type:

float

n
Type:

int

ss_order
Type:

{0, 1}

logT

log of trajectory length; used in compactification of spline coordinates

Type:

float

upper_bc_type

the boundary condition for the upper end of the spline. For ss_order == 0, this sets the derivative to zero (MSD should plateau at infinity), for ss_order == 1 we use natural boundary conditions, i.e. vanishing second derivative.

Type:

str or tuple

x_last

the theoretical bound on the x coordinate for the spline points; 2 for ss_order < 1, 1 for ss_order == 1.

Type:

{1, 2}

x_max

the maximum value for x such that the corresponding dt is still numerically representable

Type:

float

prev_fit

same as parameter previous_NPXFit_and_result.

Type:

tuple

Notes

This fit uses the same compactification scheme as SplineFit for the coordinates of the spline nodes.

Technically, the spline used to extend the MSD has n+1 nodes, since of course there has to be one at the transition from powerlaw to spline.

The vertical position of the first, and horizontal position of the last spline points are fixed. So the parameters for the spline are x0 through x(n-1) and y1 through yn.

The theoretical upper bound for the exponent of the powerlaw part is 2; at this point the covariance matrix of the process stops being positive definite. Due to numerical inaccuracies, this can start being an issue for exponents close to (but smaller than) 2 as well. Thus, the upper bound for these exponents is set to 1.99, i.e. if a fit returns 1.99 this is just the upper bound of the parameter space and not a precise estimate.

logprior(params)

Prior over the parameters

Use this function if you want to specify a prior over the parameters that will be added to the log-likelihood. Default is a flat prior, i.e. return 0.

Parameters:

params (dict) – the current parameter values; note that only the free parameters are given

Returns:

float

compactify(dt)

Compactification used for the current fit

Parameters:

dt (np.array, dtype=float) – the time lags to calculate compactification for. Should be > 0 but might include np.inf.

See also

NPXFit, decompactify_log

decompactify_log(x)

Decompactify spline points (convenience function)

Parameters:

x (np.array, dtype=float) – the compactified x-coordinates

Returns:

log_dt (np.array) – the corresponding log(dt [frames]) values

See also

compactify

params2msdm(params)

Calculate the current MSD from given parameters

Parameters:

params (dict)

Returns:

list of (msd, m)

See also

Fit.params2msdm

initial_params()

Provide an initial parameter guess

These guesses come from a powerlaw curve fit to the empirical MSD. If n > 0 we set x0 = 0.5 and insert evenly spaced spline points above that.

If previous_NPXFit_and_result was specified at initialization, we copy the powerlaw and x0 and then place evenly spaced spline points above, reproducing the existing best fit.

Returns:

params (dict)

initial_offset()

Log-likelihood offset associated with initial parameters

See Notes section of class documentation.

Returns:

float

See also

Fit

constraint_dx(params)

Make sure the spline points are properly ordered in x

We impose this constraint mainly to avoid crossing of spline points, which usually leads to the spline diverging. On top of that, conceptually this makes the solution well-defined.

Parameters:

params (np.ndarray, dtype=float) – the current fit parameters

Returns:

float – the constraint score

See also

Fit

constraint_logmsd(params)

Make sure the Spline does not diverge

Parameters:

params (np.ndarray, dtype=float) – the current fit parameters

Returns:

float – the constraint score

See also

Fit

class bayesmsd.lib.TwoLocusRouseFit(data, motion_blur_f=0)

Bases: Fit

Fit a Rouse model for two loci on a polymer at fixed separation

This class implements a fit for two loci at fixed separation, but on the same polymer. A simple model for these dynamics is given by the infinite continuous Rouse model, which gives an analytical expression for the MSD (implemented in rouse.twoLocusMSD). Here we provide an implementation to fit this MSD to data.

Parameters:

  • data (noctiluca.TaggedSet, pandas.DataFrame, list of numpy.ndarray) – the data to fit. See Fit.

  • motion_blur_f (float in [0, 1]) – fractional exposure time (i.e. exposure time as fraction of the time between frames). The resulting motion blur will be taken into account in the fit.

logprior(params)

Prior over the parameters

Use this function if you want to specify a prior over the parameters that will be added to the log-likelihood. Default is a flat prior, i.e. return 0.

Parameters:

params (dict) – the current parameter values; note that only the free parameters are given

Returns:

float

params2msdm(params)

Give an MSD function (and mean = 0) for given parameters

Parameters:

params (dict)

Returns:

list of (msd, m)

See also

rouse.twoLocusMSD

initial_params()

Initial parameters from curve fit to empirical MSD

Returns:

params (dict)

class bayesmsd.lib.DiscreteRouseFit(data, motion_blur_f=0, use_approx=False)

Bases: Fit

Fit a Rouse model for a single monomer of an infinite discrete chain (with localization error).

Parameters:

  • data (noctiluca.TaggedSet, pandas.DataFrame, list of numpy.ndarray) – the data to fit. See Fit.

  • motion_blur_f (float in [0, 1], optional) – fractional exposure time (i.e. exposure time as fraction of the time between frames). The resulting motion blur will be taken into account in the fit.

  • use_approx (bool, optional) – use a numerically more stable approximation instead of the exact expression. The approximation matches the asymptotics exactly and around the crossover stays correct to within +/-2% of the exact expression. See Notes for more details.

Notes

The analytical expression for this MSD is

\[\text{MSD}(\Delta t) = 2D\Delta t \exp\left(-\lambda\Delta t\right) \left[ I_0\left(\lambda\Delta t\right) + I_1\left(\lambda\Delta t\right) \right]\,,\]

where \(D\) is the diffusion constant of a single monomer, \(I_\nu(z)\) are Bessel functions of the first kind, \(\lambda = \frac{8D^2}{\pi\Gamma^2}\) and \(\Gamma\) is the prefactor of the large time asymptote.

Asymptotically, the above expression gives \(2D\Delta t\) at short times (the MSD of a free monomer) and \(\Gamma\sqrt{\Delta t}\) at long times (the MSD of a locus on a Rouse chain).

The exact expression above is evaluated using scipy.special.ive for the exponentially scaled Bessel functions, which gives valid results “only” for \(z < 10^9\). However, the MSD written above turns out to be represented quite well by a simple soft-min approximation between the asymptotes:

\[\text{MSD}_\text{approx.}(\Delta t) = \left[ \left(2D\Delta t\right)^{-n} + \left(\Gamma\sqrt{\Delta t}\right)^{-n} \right]^{-\frac{1}{n}}\]

with \(n = 2.5\). Clearly, this approximation is designed to match the asymptotes of the exact expression; around the crossover inbetween it stays correct to within \(\pm 2\%\). Use use_approx = True to use this approximation instead of the exact expression.

logprior(params)

Prior over the parameters

Use this function if you want to specify a prior over the parameters that will be added to the log-likelihood. Default is a flat prior, i.e. return 0.

Parameters:

params (dict) – the current parameter values; note that only the free parameters are given

Returns:

float

params2msdm(params)

Definition of MSD in terms of parameters

This is the core of the fit definition. It should give a list of tuples as required by GP.ds_logL

Parameters:

params (dict) – the current parameter values

Returns:

list of tuples (msd, m) – for the definition of msd, use of the MSDfun decorator (and potentially also imaging) is recommended.

See also

GP.ds_logL, MSDfun, imaging

initial_params()

Initial parameters from curve fit to empirical MSD

Returns:

params (dict)

class bayesmsd.lib.NPFit(data, motion_blur_f=0)

Bases: Fit

Fit a powerlaw plus noise

This is a simpler version of NPXFit, getting rid of the spline part. Note that it also uses a slightly different parametrization: (log(αΓ), α) instead of (Γ, α). This is because the former often behaves better in cases where α is not identifiable (low). In that case this new parametrization keeps Γ constant and thus identifiable.

Parameters:

  • data (noctiluca.TaggedSet, pandas.DataFrame, list of numpy.ndarray) – the data to fit. See Fit.

  • motion_blur_f (float in [0, 1], optional) – fractional exposure time (i.e. exposure time as fraction of the time between frames). The resulting motion blur will be taken into account in the fit.

logprior(params)

Prior over the parameters

Use this function if you want to specify a prior over the parameters that will be added to the log-likelihood. Default is a flat prior, i.e. return 0.

Parameters:

params (dict) – the current parameter values; note that only the free parameters are given

Returns:

float

params2msdm(params)

Definition of MSD in terms of parameters

This is the core of the fit definition. It should give a list of tuples as required by GP.ds_logL

Parameters:

params (dict) – the current parameter values

Returns:

list of tuples (msd, m) – for the definition of msd, use of the MSDfun decorator (and potentially also imaging) is recommended.

See also

GP.ds_logL, MSDfun, imaging

initial_params()

Give initial values for the parameters

You can use self.data to perform some ad hoc estimation (e.g. from the empirical MSD, using noctiluca.analysis.MSD) or just return constants.

Returns:

params (dict) – initial values for the parameters. Should have the same keys as self.parameters.

bayesmsd.parameters

class bayesmsd.parameters.Parameter(bounds=(-inf, inf), fix_to=None, linearization=None, max_linearization_moves=(10, 10))

Bases: object

Definition of one fit parameter

Parameters:

  • bounds ((lower, upper), optional) – the bounds to impose on this parameter

  • fix_to (float, str, callable, or None; optional) –

    whether and how to tie the value of this parameter to another one. Can be

    • float constant: parameter will be fixed to this value

    • name of another parameter: parameter will always be the same as the other one

    • callable: signature should be fun(params)->float, where params will be a dict with all the already determined parameter values.

    • None: keep this parameter independent.

  • linearization (Linearize.ABC, optional) – specifies a coordinate transform to a space where taking steps of 1 is (in some sense) reasonable. If unspecified, use suggest_linearization() to find a suitable transform.

  • max_linearization_moves (2-tuple of int, optional) – how far to move in the linearized space before considering the parameter unidentifiable. These are the maximum bounds in which the Profiler will explore the posterior over this parameter.

Notes

Copying a Parameter: use copy.deepcopy(). Through the linearization, a Parameter contains a self-reference; deepcopy() is smart enough to detect this and copy properly, i.e. as self-reference to the new object.

suggest_linearization()

Suggest a good linearization to use

Returns:

Linearize.ABC

class bayesmsd.parameters.Linearize

Bases: object

class ABC

Bases: object

Abstract base class for linearizations

param

the Parameter this linearization applies to. This is set by the Parameter initializer that instances of this class should be handed to.

Type:

Parameter

abstract from_linear(pe, n)

Convert from linear space to “real” parameter space

Parameters:

  • pe (float) – current point estimate

  • n (float) – coordinate in the linear space

Returns:

x (float) – coordinate in original parameter space

See also

to_linear

Notes

When subclassing, ensure that from_linear(pe, to_linear(pe, x)) == x and to_linear(pe, from_linear(pe, n)) == n for generic x and n.

from_linear should be vectorized in n, i.e. work properly on arrays.

abstract to_linear(pe, x)

Convert from “real” parameter space to linear space

Parameters:

  • pe (float) – current point estimate

  • x (float) – coordinate in original parameter space

Returns:

n (float) – coordinate in the linear space

See also

from_linear

Notes

When subclassing, ensure that from_linear(pe, to_linear(pe, x)) == x and to_linear(pe, from_linear(pe, n)) == n for generic x and n.

to_linear should be vectorized in x, i.e. work properly on arrays.

move(pe, x, n)

Move a given position x in real space by n units in the linearized space.

distance(pe, x, y)

Calculate linearized distance between x and y

mean(pe, x)

Calculate linearized mean of x (an array)

class Bounded(n_steps=10)

Bases: ABC

Linearization to use when both bounds are finite

This linearization is a simple rescaling of the “real” parameter space, such that the interval between the two bounds corresponds to n_steps units.

n = 0 corresponds to the lower bound.

Parameters:

n_steps (int) – in how many bins to subdivide the interval between the bounds

from_linear(pe, n)

Convert from linear space to “real” parameter space

Parameters:

  • pe (float) – current point estimate

  • n (float) – coordinate in the linear space

Returns:

x (float) – coordinate in original parameter space

See also

to_linear

Notes

When subclassing, ensure that from_linear(pe, to_linear(pe, x)) == x and to_linear(pe, from_linear(pe, n)) == n for generic x and n.

from_linear should be vectorized in n, i.e. work properly on arrays.

to_linear(pe, x)

Convert from “real” parameter space to linear space

Parameters:

  • pe (float) – current point estimate

  • x (float) – coordinate in original parameter space

Returns:

n (float) – coordinate in the linear space

See also

from_linear

Notes

When subclassing, ensure that from_linear(pe, to_linear(pe, x)) == x and to_linear(pe, from_linear(pe, n)) == n for generic x and n.

to_linear should be vectorized in x, i.e. work properly on arrays.

class Multiplicative(factor=2)

Bases: ABC

Log-spacing; appropriate for strictly positive/negative parameters

Each step multiplies the parameter by a constant factor. Clearly this will never cross zero, so appropriate only for strictly positive or negative parameters.

n = 0 corresponds to the point estimate.

Parameters:

factor (float) – the factor to apply at each step

from_linear(pe, n)

Convert from linear space to “real” parameter space

Parameters:

  • pe (float) – current point estimate

  • n (float) – coordinate in the linear space

Returns:

x (float) – coordinate in original parameter space

See also

to_linear

Notes

When subclassing, ensure that from_linear(pe, to_linear(pe, x)) == x and to_linear(pe, from_linear(pe, n)) == n for generic x and n.

from_linear should be vectorized in n, i.e. work properly on arrays.

to_linear(pe, x)

Convert from “real” parameter space to linear space

Parameters:

  • pe (float) – current point estimate

  • x (float) – coordinate in original parameter space

Returns:

n (float) – coordinate in the linear space

See also

from_linear

Notes

When subclassing, ensure that from_linear(pe, to_linear(pe, x)) == x and to_linear(pe, from_linear(pe, n)) == n for generic x and n.

to_linear should be vectorized in x, i.e. work properly on arrays.

class Exponential(step=1, base=2)

Bases: ABC

Take exponentially growing steps from the point estimate.

This is similar to Multiplicative in its applicability to variables that might span several orders of magnitude, but is not restricted to parameters with definite sign.

n = 0 corresponds to the point estimate; |n| = 1 corresponds to one step from the point estimate; higher n grow exponentially with given base.

Parameters:

  • step (float) – the step size at |n| = 1.

  • base (float) – the base to use for exponential growth.

from_linear(pe, n)

Convert from linear space to “real” parameter space

Parameters:

  • pe (float) – current point estimate

  • n (float) – coordinate in the linear space

Returns:

x (float) – coordinate in original parameter space

See also

to_linear

Notes

When subclassing, ensure that from_linear(pe, to_linear(pe, x)) == x and to_linear(pe, from_linear(pe, n)) == n for generic x and n.

from_linear should be vectorized in n, i.e. work properly on arrays.

to_linear(pe, x)

Convert from “real” parameter space to linear space

Parameters:

  • pe (float) – current point estimate

  • x (float) – coordinate in original parameter space

Returns:

n (float) – coordinate in the linear space

See also

from_linear

Notes

When subclassing, ensure that from_linear(pe, to_linear(pe, x)) == x and to_linear(pe, from_linear(pe, n)) == n for generic x and n.

to_linear should be vectorized in x, i.e. work properly on arrays.

bayesmsd.profiler

Implementation of the Profiler class

This is used to “explore the posterior” after running a Fit, specifically to find the credible intervals.

Example

>>> fit = <some bayesmsd.Fit subclass, c.f. bayesmsd.lib>
... profiler = bayesmsd.Profiler(fit)
... mci = profiler.find_MCI() # mci = max posterior estimate & credible
...                           # interval boundaries for each parameter

See also

bayesmsd, Profiler, Fit

class bayesmsd.profiler.Profiler(fit, profiling=True, conf=0.95, conf_tol=0.001, max_fit_runs=100, max_restarts=10, verbosity=1)

Bases: object

Exploration of the posterior after finding the point estimate

This class provides a layer on top of Fit, enabling more comprehensive exploration of the posterior after finding the (MAP) point estimate. Generally, it operates in two modes:

  • conditional posterior: wiggle each individual parameter, keeping all others fixed to the point estimate values, thus calculating conditional posterior values. In parameter space, this is a (multi-dimensional) cross along the coordinate axes.

  • profile posterior: instead of simply evaluating the posterior at each wiggle, keep the parameter fixed at the new value and optimize all others. Thus, in parameter space we are following the ridges of the posterior.

From a Bayesian point of view, beyond conditional and profile posterior, the actually interesting quantity is of course the marginal posterior. This is best obtained by sampling the posterior by MCMC (see Examples & Tutorials). Still, conditional or profile posterior often give a useful overview over the shape of the posterior. The conditional posterior is also great for setting MCMC step sizes. Profile posteriors are of course significantly more expensive than conditionals.

At the end of the day, this class moves along either profile or conditional posterior until it drops below a given cutoff, and then gives the lower bound, best value, and upper bound for each parameter. Inspired by frequentist confidence intervals, we determine the cutoff from a “confidence level” 1-α.

We follow a two-step procedure to find the lower and upper bounds: first, we move out from the point estimate with a fixed step size (in the space defined by Parameter.linearization) until the posterior drops below the cutoff. Within the thus defined bracket, we then find the actual bound to the desired accuracy by bisection. If the first step fails (i.e. the posterior does not drop below the cutoff far from the point estimate) the corresponding parameter direction is unidentifiable and the credible interval edge will be set to infinity.

Note that this class will also take care of the initial point estimate, if necessary.

The main point of entry for the user is the find_MCI method, which performs the sweeps described above. Further, run_fit might be useful if you just need the point estimate.

Parameters:

  • fit (Fit) – the Fit object to use

  • profiling (bool) – whether to run in profiling (True) or conditional (False) mode.

  • conf (float) – the “confidence level” to use for the likelihood ratio cutoff

  • conf_tol (float) – acceptable tolerance in the confidence level

  • max_fit_runs (int) – an upper bound for how often to re-run the fit (when profiling).

  • max_restarts (int) – sometimes the initial fit to find the point estimate might not converge properly, such that a better point estimate is found during the profiling runs. If that happens, the whole procedure is restarted from the new point estimate. This variable provides an upper bound on the number of these restarts. See also restart_on_better_point_estimate below.

  • verbosity ({0, 1, 2, 3}) – controls amount of messages during profiling. 0: nothing; 1: warnings only; 2: informational; 3: debugging

fit

see Parameters

Type:

Fit

min_target_from_fit

the minimization target of self.fit

Type:

Fit.MinTarget

ress

storage of all evaluated parameter points. Each entry corresponds to a parameter dimension and contains a list of dicts that were obtained while sweeping that parameter. The individual entries are dicts like the one returned by Fit.run, with fields 'params' and 'logL'.

Type:

dict

point_estimate

the MAP point estimate, a dict like the other points in ress.

Type:

dict

conf, conf_tol

see Parameters

Type:

float

LR_target

the target decline in the posterior value we’re searching for

Type:

float

LR_interval

acceptable interval corresponding to conf_tol

Type:

[float, float]

cur_param

the name of the parameter currently being sweeped.

Type:

str

max_fit_runs

see Parameters

Type:

int

run_count

counts the runs executed so far

Type:

int

max_restarts_per_parameters

same as max_restarts parameter.

Type:

int

verbosity

see Parameters

Type:

int

restart_on_better_point_estimate

whether to restart upon finding a better point estimate. Might make sense to disable if your posterior is rugged, to avoid restarting all the time.

Type:

bool

bar

the progress bar showing successive fit evaluations

Type:

tqdm progress bar

See also

run_fit, find_MCI

vprint(verbosity, *args, **kwargs)

Prints only if self.verbosity >= verbosity.

property profiling
exception FoundBetterPointEstimate

Bases: Exception

restart_if_better_pe_found()

Decorator taking care of restarts

We use the Profiler.FoundBetterPointEstimate exception to handle restarts. So a function decorated with this decorator can just raise this exception and will then be restarted properly.

static likelihood_significantly_greater(res1, res2)

Helper function

The threshold is 1e-3. Note that this is an asymmetric operation, i.e. there is a regime where neither res1 significantly greater res2 nor the other way round.

Parameters:

  • res1 (dict) – like the output of Fit.run, and the entries of self.ress.

  • res2 (dict) – like the output of Fit.run, and the entries of self.ress.

Returns:

bool

property best_estimate

The best current estimate

This should usually be the point estimate, but we might find a better one along the way.

check_point_estimate_against(res)

Check whether res is better than current point estimate

Parameters:

res (dict) – the evaluated parameter point to check. A dict like the ones in self.ress.

Raises:

Profiler.FoundBetterPointEstimate –

run_fit(is_new_point_estimate=True, **fit_kw)

Execute one fit run

This is used to find the initial point estimate, as well as subsequent evaluations of the profile posterior.

Parameters:

  • is_new_point_estimate (bool) – whether we are looking for a new point estimate or just evaluating a profile point. This just affects where the result is stored internally

  • fit_kw (keyword arguments) – additional parameters for Fit.run

See also

find_MCI

find_closest_res(val, direction=None)

Find the closest previously evaluated point

Parameters:

  • val (float) – new value

  • direction ({None, -1, 1}) – search only for existing values that are greater (1) or smaller (-1) than the specified one.

Returns:

dict – appropriate point from self.ress

See also

profile_likelihood

profile_likelihood(value, init_from='closest')

Evaluate profile (or conditional) likelihood / posterior

Parameters:

  • value (float) – value of the current parameter of interest (self.cur_param)

  • init_from (dict or 'closest') – from where to start the optimization (only relevant if self.profiling). Set to ‘closest’ to use self.find_closest_res.

Returns:

dict – like in self.ress (and also stored there)

find_bracket_point(direction)

First step of the procedure (“bracketing”)

In this first step we simply push out from the point estimate to try and establish a bracket containing the exact boundary point. See Profiler.bracket_strategy. This function implements that push, for one parameter in one direction.

Parameters:

direction ({-1, 1}) – which direction to go in

Returns:

x, pL (float) – the found bracket end point and associated posterior. If the chosen direction turns out to be unidentifiable, x is set to the corresponding bound (might be inf), and pL = np.inf.

See also

initial_bracket_points

initial_bracket_points()

Execute bracketing in both directions

See also

find_bracket_point, find_MCI

solve_bisection(bracket, bracket_pL)

Find exact root by recursive bisection

Parameters:

  • lin_bracket ([float, float]) – the bracket of parameter values containing the root.

  • bracket_pL ([float, float]) – the posterior values associated with the bracket points

Returns:

float – the root within the bracket, to within the given precision (c.f self.conf_tolerance)

Notes

All points evaluated along the way are stored in self.ress

find_single_MCI(param_name)

Run the whole profiling process for one parameter

Parameters:

param_name (str) – name of the parameter to sweep

Returns:

  • m (float) – the point estimate value for this parameter

  • ci (np.array([float, float])) – the bounds where the posterior has dropped below the maximum value the specified amount (c.f. self.conf)

See also

find_MCI

find_MCI(parameters=<built-in function all>, show_progress=False)

Perform the sweep for multiple parameters

Parameters:

parameters (str or list of str, optional) – the names of the parameters to sweep. Defaults to “all independent parameters from the fit”

Returns:

dict – for each parameter a tuple (m, ci), where m is the point estimate and ci a numpy array with lower and upper bound; c.f. find_single_MCI

See also

run_fit, find_single_MCI