Bout analysis

Tools and classes for the identification of behavioural bouts

A histogram of log-transformed frequencies of x with a chosen bin width and upper limit forms the basis for models. Histogram bins following empty ones have their frequencies averaged over the number of previous empty bins plus one. Models attempt to discern the number of random Poisson processes, and their parameters, generating the underlying distribution of log-transformed frequencies.

The abstract class Bouts provides basic methods.

Abstract class & methods summary

Bouts(x, bw[, method])	Abstract base class for models of log-transformed frequencies
Bouts.init_pars(x_break[, plot, ax])	Find starting values for mixtures of random Poisson processes
Bouts.fit(start, **kwargs)	Fit Poisson mixture model to log frequencies
Bouts.bec(coefs)	Calculate bout ending criteria from model coefficients
Bouts.plot_fit(coefs[, ax])	Plot log frequency histogram and fitted model

Nonlinear least squares models

Currently, the model describing the histogram as it is built is implemented in the BoutsNLS class. For the case of a mixture of two Poisson processes, this class would set up the model:

(1)\[y = log[N_f \lambda_f e^{-\lambda_f t} + N_s \lambda_s e^{-\lambda_s t}]\]

where \(N_f\) and \(N_s\) are the number of events belonging to process \(f\) and \(s\), respectively; and \(\lambda_f\) and \(\lambda_s\) are the probabilities of an event occurring in each process. Mixtures of more processes can also be added to the model.

The bout-ending criterion (BEC) corresponding to equation (1) is:

(2)\[BEC = \frac{1}{\lambda_f - \lambda_s} log \frac{N_f \lambda_f}{N_s \lambda_s}\]

Note that there is one BEC per transition between Poisson processes.

The methods of this subclass are provided by the abstract super class Bouts, and defining those listed below.

Methods summary

BoutsNLS.plot_ecdf(coefs[, ax])

Plot observed and modelled empirical cumulative frequencies

Maximum likelihood models

This is the preferred approach to modelling mixtures of random Poisson processes, as it does not rely on the subjective construction of a histogram. The histogram is only used to generate reasonable starting values, but the underlying paramters of the model are obtained via maximum likelihood, so it is more robust.

For the case of a mixture of two processes, as above, the log likelihood of all the \(N_t\) in a mixture can be expressed as:

(3)\[log\ L_2 = \sum_{i=1}^{N_t} log[p \lambda_f e^{-\lambda_f t_i} + (1-p) \lambda_s e^{-\lambda_s t_i}]\]

where \(p\) is a mixing parameter indicating the proportion of fast to slow process events in the sampled population.

The BEC in this case can be estimated as:

(4)\[BEC = \frac{1}{\lambda_f - \lambda_s} log \frac{p\lambda_f}{(1-p)\lambda_s}\]

The subclass BoutsMLE offers the framework for these models.

Methods summary

BoutsMLE.negMLEll(params, x[, istransformed])	Log likelihood function of parameters given observed data
BoutsMLE.fit(start[, fit1_opts, fit2_opts])	Maximum likelihood estimation of log frequencies
BoutsMLE.bec(fit)	Calculate bout ending criteria from model coefficients
BoutsMLE.plot_fit(fit[, ax])	Plot log frequency histogram and fitted model
BoutsMLE.plot_ecdf(fit[, ax])	Plot observed and modelled empirical cumulative frequencies

class bouts.Bouts(x, bw, method='standard')

Bases: object

Abstract base class for models of log-transformed frequencies

This is a base class for other classes to build on, and do the model fitting. Bouts is an abstract base class to set up bout identification procedures. Subclasses must implement fit and bec methods, or re-use the default NLS methods in Bouts.

x

1D array with input data.

Type:

array_like

method

Method used for calculating the histogram.

Type:

str

lnfreq

DataFrame with the centers of histogram bins, and corresponding log-frequencies of x.

Type:

pandas.DataFrame

abstract bec(coefs)

Calculate bout ending criteria from model coefficients

Implementing default as from NLS method.

Parameters:

coefs (pandas.DataFrame) – DataFrame with model coefficients in columns, and indexed by parameter names “a” and “lambda”.

Returns:

out – 1-D array with BECs implied by coefs. Length is coefs.shape[1]

Return type:

numpy.ndarray, shape (n,)

abstract fit(start, **kwargs)

Fit Poisson mixture model to log frequencies

Default is non-linear least squares method.

Parameters:

• start (pandas.DataFrame) – DataFrame with coefficients for each process in columns.

• **kwargs – Optional keyword arguments passed to scipy.optimize.curve_fit().

Returns:

• coefs (pandas.DataFrame) – Coefficients of the model.

• pcov (2D array) – Covariance of coefs.

init_pars(x_break, plot=True, ax=None, **kwargs)

Find starting values for mixtures of random Poisson processes

Starting values are calculated using the “broken stick” method.

Parameters:

• x_break (array_like) – One- or two-element array with values determining the break(s) for broken stick model, such that \(x < x\_break_0\) is first process, \(x >= x\_break_1\) & \(x < x\_break_2\) is second process, and \(x >= x\_break_2\) is third one.

• plot (bool, optional) – Whether to plot the broken stick model.

• ax (matplotlib.Axes, optional) – An Axes instance to use as target. Default is to create one.

• **kwargs (optional keyword arguments) – Passed to plotting function. Ignored currently.

Returns:

out – DataFrame with coefficients for each process.

Return type:

pandas.DataFrame

plot_fit(coefs, ax=None)

Plot log frequency histogram and fitted model

Parameters:

• coefs (pandas.DataFrame) – DataFrame with model coefficients in columns, and indexed by parameter names “a” and “lambda”.

• ax (matplotlib.Axes instance) – An Axes instance to use as target.

Returns:

ax

Return type:

matplotlib.Axes

class bouts.BoutsMLE(x, bw, method='standard')

Bases: Bouts

Maximum Likelihood estimation for models of Poisson process mixtures

Methods for modelling log-frequency data as a mixture of Poisson processes via maximum likelihood estimation [2], [3]. Mixtures of two or three Poisson processes are supported.

Even in these relatively simple cases, it is very important to provide good starting values for the parameters.

One useful strategy to get good starting parameter values is to proceed in 4 steps. First, fit a broken stick model to the log frequencies of binned data (see init_pars()), to obtain estimates of 4 parameters in a 2-process model [1], or 6 in a 3-process model. Second, calculate parameter(s) \(p\) from the \(\alpha\) parameters obtained by fitting the broken stick model, to get tentative initial values as in [2]. Third, obtain MLE estimates for these parameters, but using a reparameterized version of the -log L2 function. Lastly, obtain the final MLE estimates for the three parameters by using the estimates from step 3, un-transformed back to their original scales, maximizing the original parameterization of the -log L2 function.

init_pars() can be used to perform step 1. Calculation of the mixing parameters \(p\) in step 2 is trivial from these estimates. Method negMLEll() calculates the negative log-likelihood for a reparameterized version of the -log L2 function given by [1], so can be used for step 3. This uses a logit transformation of the mixing parameter \(p\), and log transformations for density parameters \(\lambda\). Method negMLEll() is used again to compute the -log L2 function corresponding to the un-transformed model for step 4.

The fit() method performs the main job of maximizing the -log L2 functions, and is essentially a wrapper around minimize(). It only takes the -log L2 function, a DataFrame of starting values, and the variable to be modelled, all of which are passed to minimize() for optimization. Additionally, any other arguments are also passed to minimize(), hence great control is provided for fitting any of the -log L2 functions.

In practice, step 3 does not pose major problems using the reparameterized -log L2 function, but it might be useful to use method ‘L-BFGS-B’ with appropriate lower and upper bounds. Step 4 can be a bit more problematic, because the parameters are usually on very different scales and there can be multiple minima. Therefore, it is almost always the rule to use method ‘L-BFGS-B’, again bounding the parameter search, as well as other settings for controlling the optimization.

References

[2] (1,2)

Langton, S.; Collett, D. and Sibly, R. (1995) Splitting behaviour into bouts; a maximum likelihood approach. Behaviour 132, 9-10.

[3]

Luque, S.P. and Guinet, C. (2007) A maximum likelihood approach for identifying dive bouts improves accuracy, precision, and objectivity. Behaviour, 144, 1315-1332.

Examples

See Simulating bouts for a detailed example.

bec(fit)

Calculate bout ending criteria from model coefficients

Parameters:

fit (scipy.optimize.OptimizeResult) – Object with the optimization result, having a x attribute with coefficients of the solution.

Returns:

out

Return type:

numpy.ndarray

Notes

Current implementation is for a two-process mixture, hence an array of a single float is returned.

fit(start, fit1_opts=None, fit2_opts=None)

Maximum likelihood estimation of log frequencies

Parameters:

• start (pandas.DataFrame) – DataFrame with starting values for coefficients of each process in columns. These can come from the “broken stick” method as in Bouts.init_pars(), and will be transformed to minimize the first log likelihood function.

• fit1_opts (dict) – Dictionaries with keywords to be passed to scipy.optimize.minimize(), for the first and second fits.

• fit2_opts (dict) – Dictionaries with keywords to be passed to scipy.optimize.minimize(), for the first and second fits.

Returns:

fit1, fit2 – Objects with the optimization result from the first and second fit, having a x attribute with coefficients of the solution.

Return type:

scipy.optimize.OptimizeResult

Notes

Current implementation handles mixtures of two Poisson processes.

negMLEll(params, x, istransformed=True)

Log likelihood function of parameters given observed data

Parameters:

• params (array_like) – 1-D array with parameters to fit. Currently must be either 3-length, with mixing parameter \(p\), density parameter \(\lambda_f\) and \(\lambda_s\), in that order, or 5-length, with \(p_f\), \(p_fs\), \(\lambda_f\), \(\lambda_m\), \(\lambda_s\).

• x (array_like) – Independent data array described by model with parameters params.

• istransformed (bool) – Whether params are transformed and need to be un-transformed to calculate the likelihood.

Return type:

out

plot_ecdf(fit, ax=None, **kwargs)

Plot observed and modelled empirical cumulative frequencies

Parameters:

• fit (scipy.optimize.OptimizeResult) – Object with the optimization result, having a x attribute with coefficients of the solution.

• ax (matplotlib.axes.Axes instance) – An Axes instance to use as target.

• **kwargs – Optional keyword arguments passed to matplotlib.pyplot.gca().

Returns:

Axes instances.

Return type:

ax

plot_fit(fit, ax=None)

Plot log frequency histogram and fitted model

Parameters:

• fit (scipy.optimize.OptimizeResult) – Object with the optimization result, having a x attribute with coefficients of the solution.

• ax (matplotlib.axes.Axes instance) – An Axes instance to use as target.

Returns:

Axes instances.

Return type:

ax

class bouts.BoutsNLS(x, bw, method='standard')

Bases: Bouts

Nonlinear Least Squares fitting for models of Poisson process mixtures

Methods for modelling log-frequency data as a mixture of Poisson processes via nonlinear least squares [1].

References

[1] (1,2,3)

Sibly, R.; Nott, H. and Fletcher, D. (1990) Splitting behaviour into bouts Animal Behaviour 39, 63-69.

Examples

Draw 1000 samples from a mixture where the first process occurs with \(p < 0.7\) and the second process occurs with the remaining probability.

>>> from skdiveMove.tests import random_mixexp
>>> rng = np.random.default_rng(123)
>>> x2 = random_mixexp(1000, p=0.7, lda=np.array([0.05, 0.005]),
...                    rng=rng)
>>> xbouts2 = BoutsNLS(x2, bw=5)
>>> init_pars = xbouts2.init_pars([80], plot=False)

Fit the model and retrieve coefficients:

>>> coefs, pcov = xbouts2.fit(init_pars)
>>> print(np.round(coefs, 4))
        (2.519, 80.0]  (80.0, 1297.52]
a           3648.8547        1103.4423
lambda         0.0388           0.0032

Calculate bout-ending criterion (returns array):

>>> print(np.round(xbouts2.bec(coefs), 4))
[103.8648]

Plot observed and predicted data:

>>> xbouts2.plot_fit(coefs)  
<Axes: ...>

Plot ECDF:

>>> xbouts2.plot_ecdf(coefs)  
<Axes: ...>

bec(coefs)

Calculate bout ending criteria from model coefficients

The metaclass bouts.Bouts implements this method.

Parameters:

coefs (pandas.DataFrame) – DataFrame with model coefficients in columns.

Returns:

out – 1-D array with BECs implied by coefs. Length is coefs.shape[1]

Return type:

numpy.ndarray, shape (n,)

fit(start, **kwargs)

Fit non-linear least squares model to log frequencies

The metaclass bouts.Bouts implements this method.

Parameters:

• start (pandas.DataFrame) – DataFrame with coefficients for each process in columns.

• **kwargs – Optional keyword arguments passed to scipy.optimize.curve_fit().

Returns:

• coefs (pandas.DataFrame) – Coefficients of the model.

• pcov (2D array) – Covariance of coefs.

plot_ecdf(coefs, ax=None, **kwargs)

Plot observed and modelled empirical cumulative frequencies

Parameters:

• coefs (pandas.DataFrame) – DataFrame with model coefficients in columns.

• ax (matplotlib.axes.Axes instance) – An Axes instance to use as target.

• **kwargs – Optional keyword arguments passed to matplotlib.pyplot.gca().

Returns:

Axes instances.

Return type:

ax

bouts.label_bouts(x, bec, as_diff=False)

Classify data into bouts based on bout ending criteria

Parameters:

• x (pandas.Series) – Series with data to classify according to bec.

• bec (array_like) – Array with bout-ending criteria. It is assumed to be sorted.

• as_diff (bool, optional) – Whether to apply diff on x so it matches bec’s scale.

Returns:

out – Integer array with the same shape as x.

Return type:

numpy.ndarray

bouts.random_mixexp(n, p, lda, rng=None)

Generate samples from mixture of exponential distributions

Simulate a mixture of two or three random exponential distributions. This uses a special definition for the probabilities (\(p_i\)). In the two-process case, \(p\) represents the proportion of “fast” to “slow” events in the mixture. In the three-process case, \(p_0\) represents the proportion of “fast” to “slow” events, and \(p_1\) represents the proportion of “slow” to “slow” and “very slow” events.

Parameters:

• n (int) – Output sample size.

• p (float or array_like) – Probabilities for processes in the output mixture sample.

• lda (array_like) – array_like with \(\lambda\) (scale) for each process.

• rng (Generator) – Random number generator object. If not provided, a default one is used.

Return type:

ndarray

Examples

Draw 1000 samples from a mixture where the first process occurs with \(p < 0.7\) and the second process occurs with the remaining probability.

>>> rng = np.random.default_rng(123)
>>> random_mixexp(1000, p=0.7, lda=np.array([0.05, 0.005]),
...               rng=rng)  
array([7.27468930e+00, 2.69331091e+00, 5.29542594e+00, 1.02697947e+01,
       ...])