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.
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) – the effective short time scaling exponent.

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.

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.

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
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 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.

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:

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)

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.

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

run(init_from=None, optimization_steps=('simplex',), maxfev=None, fix_values=None, 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.
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).

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.

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 == 0 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, 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, 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 == 0, 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 == 0, 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

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 == 0 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 == 0, 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.

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.

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)

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

profiling

see Parameters

Type:: bool

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

API reference

bayesmsd.deco

bayesmsd.fit

bayesmsd.gp

bayesmsd.lib

bayesmsd.parameters

bayesmsd.profiler

bayesmsd.deco 

bayesmsd.fit 

bayesmsd.gp 

bayesmsd.lib 

bayesmsd.parameters 

bayesmsd.profiler 