birdepy.estimate()

birdepy.estimate(t_data, p_data, p0, p_bounds, framework='dnm', model='Verhulst', scheme='discrete', con=(), known_p=(), idx_known_p=(), se_type='asymptotic', seed=None, ci_plot=False, export=False, display=False, **options)

Parameter estimation for (continuously or discretely observed) continuous-time birth-and-death processes.

Parameters:

  • t_data (list) – Observation times of birth-and-death process. If one trajectory is observed, then this is a list. If multiple trajectories are observed, then this is a list of lists where each list corresponds to a trajectory.

  • p_data (list) – Observed populations of birth-and-death process at times in argument t_data. If one trajectory is observed, then this is a list. If multiple trajectories are observed, then this is a list of lists where each list corresponds to a trajectory.

  • p0 (array_like) – Initial parameter guess. Array of real elements of size (n,), where n is the number of unknown parameters.

  • p_bounds (list) – Bounds on parameters. Should be specified as a sequence of (min, max) pairs for each unknown parameter. See here.

  • framework (string, optional) –

    Parameter estimation framework. Should be one of:

    • ’abc’ (see here)

    • ’dnm’ (default) (see here)

    • ’em’ (see here)

    • ’lse’ (see here)

    Additional kwargs are available for each framework.

  • model (string, optional) –

    Model specifying birth and death rates of process (see here). Should be one of:

    • ’Verhulst’ (default)

    • ’Ricker’

    • ’Hassell’

    • ’MS-S’

    • ’Moran’

    • ’pure-birth’

    • ’pure-death’

    • ’Poisson’

    • ’linear’

    • ’linear-migration’

    • ’M/M/1’

    • ’M/M/inf’

    • ’loss-system’

    • ’custom’

    If set to ‘custom’, then kwargs b_rate and d_rate must also be specified. See here for more information.

  • scheme (string, optional) –

    Observation scheme. Should be one of:

    • ’discrete’ (default)

    • ’continuous’ (see here)

    If set to ‘continuous’, then it is assumed that the population is observed continuously with jumps occuring at times in t_data into corresponding states in p_data.

  • con ({Constraint, dict} or List of {Constraint, dict}, optional) – Constraints definition for parameters (only for kwarg opt_method equal to ‘differential-evolution’, ‘COBYLA’, ‘SLSQP’ or ‘trust-constr’). See here for more information.

  • known_p (array_like, optional) – List of known parameter values. For built in models these must be in their canonical order as given (here). If this argument is given, then argument idx_known_p must also be specified. See here for more information.

  • idx_known_p (array_like, optional) – List of indices of known parameters (as given in argument ‘known_p’). For built in models indices must correspond to canonical order as given here. If this argument is given, then argument known_p must also be specified. See here for more information.

  • opt_method (str or callable, optional) –

    Type of solver. Should be one of the options available to scipy.optimize.minimize() (see here for more information) or ‘differential-evolution’ (see here for more information).

  • se_type (string, optional) – Should be one of: ‘none’ (default), ‘simulated’, or ‘asymptotic’. See here for more information.

  • display (bool, optional) – If True, then progress updates are provided for some methods.

  • ci_plot (bool, optional) – Enables confidence region plotting for 2d parameter estimates.

  • export (str, optional) – File name for export of confidence region figure to a LaTeX file.

Returns:

res – The estimation output represented as an EstimationOutput() object. Important attributes are: p the parameter estimate, se the standard error estimate, cov the estimated covariance of the assumed distribution of the parameter estimate, val is the log-likelihood for ‘framework’ dnm’ and ‘em’, squared error for ‘framework’ `lse’, `capacity is the estimated carrying capacity.

Return type:

EstimationOutput

Examples

Example 1: Simulate a discretely observed sample path and estimate the parameters using the alternative frameworks. First simulate some sample paths of a Ricker model using birdepy.simulate.discrete():

import birdepy as bd
t_data = list(range(100))
p_data = bd.simulate.discrete([0.75, 0.25, 0.02, 1], 'Ricker', 10, t_data,
                              survival=True, seed=2021)

Then, using the simulated data, estimate the parameters using ‘dnm’, ‘lse’ and ‘em’ as the argument of framework:

est_dnm = bd.estimate(t_data, p_data, [0.5, 0.5, 0.05], [[0,1], [0,1], [0, 0.1]],
                      framework='dnm', model='Ricker', idx_known_p=[3], known_p=[1])
est_em = bd.estimate(t_data, p_data, [1, 0.5, 0.05], [[0,1], [0,1], [1e-6,0.1]],
                      framework='em', model='Ricker', idx_known_p=[3], known_p=[1])
est_abc = bd.estimate(t_data, p_data, [0.5, 0.5, 0.05], [[0,1], [0,1], [0, 0.1]],
                      framework='abc', model='Ricker', idx_known_p=[3], known_p=[1])
est_lse = bd.estimate(t_data, p_data, [0.5, 0.5, 0.05], [[0,1], [0,1], [0, 0.1]],
                      framework='lse', model='Ricker', idx_known_p=[3], known_p=[1], se_type='simulated')
print(f'dnm estimate: {est_dnm.p}, dnm standard errors: {est_dnm.se}')
print(f'lse estimate: {est_lse.p}, lse standard errors: {est_lse.se}')
print(f'abc estimate: {est_abc.p}, abc standard errors: {est_abc.se}')
print(f'em estimate: {est_em.p}, em standard errors: {est_em.se}')

Outputs:

dnm estimate: [0.7477212189824904, 0.2150484536334751, 0.022745124483227304] , dnm standard errors: [0.16904225 0.03443199 0.00433567]
em estimate: [0.7375802511179848, 0.19413965548145604, 0.024402343633644553] , em standard errors: [0.15742852 0.02917437 0.00429763]
abc estimate: [0.6318632898413052, 0.02074882329749562, 0.06580340596326038], abc standard errors: [0.22865334, 0.0148124, 0.0129306]
lse estimate: [0.7941741586214265, 0.2767698457541133, 0.01935636627568731] , lse standard errors: [0.1568291  0.19470746 0.01243208]

Alternatively, we may be interested in continuously observed data.

Example 2: Simulate a continuous sample path and estimate the parameters.

Simulate some synthetic data:

t_data, p_data = bd.simulate.continuous([0.75, 0.25, 0.02, 1], 'Ricker', 10,
                                        100, survival=True, seed=2021)

Estimate:

est = bd.estimate(t_data, p_data, [0.5, 0.5, 0.05], [[0,1], [0,1], [0, 0.1]],
                  model='Ricker', idx_known_p=[3], known_p=[1], scheme='continuous')
print(est.p)

Outputs:

[0.7603171062895576, 0.2514810854871476, 0.020294342655751033]

Notes

For practitioners aiming to fit a specific model to their data, the default methods for estimation should generally suffice to meet their needs. However, it is advisable for practitioners to evaluate the estimation results using their expert judgment. If the outcomes do not align with existing theory or other established facts, it is recommended that they explore alternative options. For example, setting the opt_method parameter to ‘differential-evolution’ may provide more accurate results at the cost of longer computation times.

Estimation algorithms and models are also described in [1]. If you use this function for published work, then please cite this paper.

For a text book treatment on the theory of birth-and-death processes see [2].

If model is ‘Poisson’, then the true dnm is immediately returned (i.e., the total number of events during the observation periods divided by the total observation time).

See also

birdepy.estimate() birdepy.probability() birdepy.forecast()

birdepy.simulate.discrete() birdepy.simulate.continuous()

birdepy.gpu_functions.probability() birdepy.gpu_functions.discrete()

References

[View birdepy.estimate() source code]

birdepy.estimate(framework=’abc’)

birdepy.interface_abc.discrete_est_abc(data, eps_abc, distance, stat, k, gam, max_its, max_q, eps_change, method, b_rate, d_rate, idx_known_p, known_p, p_bounds, con, tau, seed, display)

Parameter estimation for discretely observed continuous-time birth-and-death processes using approximate Bayesian computation. See here for more information.

To use this function call birdepy.estimate() with framework set to ‘abc’:

bd.estimate(t_data, p_data, p0, p_bounds, framework='abc', eps_abc='dynamic', k=100,
            max_its=3, max_q=0.99, eps_change=5, gam=5, method='gwa', tau=None, seed=None,
            distance=None, stat='mean', display=False)

The parameters associated with this framework (listed below) can be accessed using kwargs in birdepy.estimate(). See documentation of birdepy.estimate() for the main arguments.

Parameters:

  • eps_abc (list, str, optional) – Threshold distance betweeen simulated data and observed data for accepting parameter proposals. If set to ‘dynamic’ (default), then the procedure described in [3] is used. Otherwise eps_abc must be a list which specifies epsilon for each iteration.

  • k (int, optional) – Number of successful parameter samples used to obtain estimate.

  • max_its (int, optional) – Maximum number of iterations of algorithm.

  • max_q (scalar, optional) – Tolerance threshold for stopping algorithm (see Equation 2.5 in [3]). Is only checked after at least two iterations have occurred.

  • eps_change (scalar, optional) – An iteration is only performed if the percentage decrease in ‘eps_abc’ compared to the previous iteration is greater than this value.

  • gam (int, optional) – If eps_abc is set to ‘dynamic’, then k*gam samples are initially sampled and the distance between the data and the k-th largest of these samples is used as the first value of epsilon

  • method (string, optional) –

    Simulation algorithm used to generate samples (see here). Should be one of:

    • ’exact’

    • ’ea’

    • ’ma’

    • ’gwa’ (default)

  • tau (scalar, optional) – Time between samples for the approximation methods ‘ea’, ‘ma’ and ‘gwa’. Has default value min(x/10, 0.1) where ‘x’ is the smallest inter-observation time.

  • seed (int, Generator, optional) – If seed is not specified the random numbers are generated according to np.random.default_rng(). If seed is an ‘int;, random numbers are generated according to np.random.default_rng(seed)(). If seed is a ‘Generator’, then that object is used. See here for more information.

  • distance (callable, optional) – Computes the distance between simulated data and observed data. Default value is scipy.spatial.distance.euclidean().

  • stat (string, optional) – Determines which statistic is used to summarize the posterior distribution. Should be one of: ‘mean’ or ‘median’.

Examples

Simulate a sample path and estimate the parameters using ABC:

import birdepy as bd
t_data = list(range(100))
p_data = bd.simulate.discrete([0.75, 0.25, 0.02, 1], 'Ricker', 10, t_data,
                              survival=True, seed=2021)

Assume that the death rate and population size are known, then estimate the rate of spread:

est = bd.estimate(t_data, p_data, [0.5], [[0,1]], framework='abc',
                  model='Ricker', idx_known_p=[1, 2, 3],
                  known_p=[0.25, 0.02, 1], display=True, seed=2021)
print(f"abc estimate is {est.p}, with standard errors {est.se},
      computed in {est.compute_time} seconds.")

Outputs:

abc estimate is [0.72434934654868] , with standard errors [0.0520378066100896] computed in  122.48563146591187 seconds.

Notes

Estimation algorithms and models are also described in [1]. If you use this function for published work, then please cite this paper.

For a text book treatment on the theory of birth-and-death processes see [2].

See also

birdepy.estimate() birdepy.probability() birdepy.forecast()

birdepy.simulate.discrete() birdepy.simulate.continuous()

References

[View birdepy.estimate(framework=’abc’) source code]

birdepy.estimate(framework=’dnm’)

birdepy.interface_dnm.discrete_est_dnm(data, likelihood, model, z_trunc, idx_known_p, known_p, p0, p_bounds, con, opt_method, options)

Parameter estimation for discretely observed continuous-time birth-and-death processes using direct numerical maximization of approximate likelihood functions. See here for more information.

To use this function call birdepy.estimate() with framework set to ‘dnm’:

bd.estimate(t_data, p_data, p0, p_bounds, framework='dnm', likelihood='expm', z_trunc=())

The parameters associated with this framework (listed below) can be accessed using kwargs in birdepy.estimate(). See documentation of birdepy.estimate() for the main arguments.

Parameters:

  • likelihood (string, optional) –

    Likelihood approximation method. Should be one of (alphabetical ordering):

    • ’da’ (see here)

    • ’Erlang’ (see here)

    • ’expm’ (default) (see here)

    • ’gwa’ (see here)

    • ’gwasa’ (see here)

    • ’ilt’ (see here)

    • ’oua’ (see here)

    • ’uniform’ (see here)

    The links point to the documentation of the relevant method in birdepy.probability(). The arguments associated with each of these methods may be used as a kwarg in birdepy.estimate() when likelihood is set to use the method.

  • z_trunc (array_like, optional) – Truncation thresholds, i.e., minimum and maximum states of process considered. Array of real elements of size (2,) by default z_trunc=[z_min, z_max] where z_min=max(0, min(p_data) - 2*z_range) and z_max=max(p_data) + 2*z_range where z_range=max(p_data)-min(p_data). Only applicable to likelihood methods ‘Erlang’, ‘expm’ and ‘uniform’.

Examples

Simulate a sample path and estimate the parameters using the various likelihood approximation methods:

import birdepy as bd
t_data = list(range(100))
p_data = bd.simulate.discrete([0.75, 0.25, 0.02, 1], 'Ricker', 10, t_data,
                              survival=True, seed=2021)
for likelihood in ['da', 'Erlang', 'expm', 'gwa', 'gwasa', 'ilt', 'oua', 'uniform']:
    likelihood = 'gwasa'
    est = bd.estimate(t_data, p_data, [0.5, 0.5, 0.05], [[1e-6,1], [1e-6,1], [1e-6, 0.1]],
                     framework='dnm', model='Ricker', idx_known_p=[3], known_p=[1],
                      likelihood=likelihood)
    print('dnm estimate using', likelihood, 'is', est.p, ', with standard errors',
          est.se, 'computed in ', est.compute_time, 'seconds.')

Outputs:

dnm estimate using da is [0.742967758293063, 0.21582082476775125, 0.022598554340938885] , with standard errors [0.16633436 0.03430501 0.00436483] computed in  12.29400086402893 seconds.
dnm estimate using Erlang is [0.7478158070214841, 0.2150341741826537, 0.022748822721356705] , with standard errors [0.16991822 0.0345762  0.00435054] computed in  0.9690003395080566 seconds.
dnm estimate using expm is [0.7477255598292176, 0.2150476994316206, 0.022745305129350565] , with standard errors [0.16904095 0.03443197 0.00433563] computed in  1.6919987201690674 seconds.
dnm estimate using gwa is [0.6600230500097711, 0.16728663936008945, 0.02512248420514078] , with standard errors [0.14248815 0.02447161 0.00488879] computed in  37.52255415916443 seconds.
dnm estimate using gwasa is [0.6604981297820195, 0.16924607541398484, 0.02492054535741541] , with standard errors [0.14244908 0.02485465 0.00488222] computed in  0.8699958324432373 seconds.
dnm estimate using ilt is [0.7466254648849691, 0.21415145383850764, 0.022794996238547492] , with standard errors [0.10187377 0.03137803        nan] computed in  1185.0924031734467 seconds.
dnm estimate using oua is [0.5000083585920406, 0.5, 0.05] , with standard errors [       nan        nan 0.01961143] computed in  3.466001272201538 seconds.
dnm estimate using uniform is [0.7477293759434092, 0.215047068344254, 0.022745437226772615] , with standard errors [0.16900378 0.03443071 0.00433504] computed in  3.275972366333008 seconds.

The constraint con={'type': 'ineq', 'fun': lambda p: p[0]-p[1]} ensures that p[0] > p[1] (i.e., rate of spread greater than recovery rate).

How well methods perform varies from case to case. In this instance most methods perform well, while some throw errors but return useful output regardless, and some fail altogether.

Notes

Estimation algorithms and models are also described in [1]. If you use this function for published work, then please cite this paper.

For a text book treatment on the theory of birth-and-death processes see [2].

See also

birdepy.estimate() birdepy.probability() birdepy.forecast()

birdepy.simulate.discrete() birdepy.simulate.continuous()

References

[View birdepy.estimate(framework=’dnm’) source code]

birdepy.estimate(framework=’em’)

birdepy.interface_em.discrete_est_em(data, p0, technique, accelerator, likelihood, p_bounds, con, known_p, idx_known_p, model, b_rate, d_rate, z_trunc, max_it, i_tol, j_tol, h_tol, display, opt_method, options)

Parameter estimation for discretely observed continuous-time birth-and-death processes using the expectation maximization algorithm. See here for more information.

To use this function call birdepy.estimate() with framework set to ‘em’:

bd.estimate(t_data, p_data, p0, p_bounds, framework='em', technique='expm',
                 accelerator='none', likelihood='expm', laplace_method='cme-talbot',
                 lentz_eps=1e-6, max_it=25, i_tol=1e-2, j_tol=1e-1, h_tol=1e-2,
                 z_trunc=())

The parameters associated with this framework (listed below) can be accessed using kwargs in birdepy.estimate(). See documentation of birdepy.estimate() for the main arguments.

Parameters:

  • technique (string, optional) –

    Expectation step technique. Should be one of (alphabetical ordering):

    • ’expm’ (default)

    • ’ilt’

    • ’num’

    See here for more information.

  • accelerator (string, optional) –

    EM accelerator method. Should be one of (alphabetical ordering):

    • ’cg’ (see [3])

    • ’none’ (see [4])

    • ’Lange’ (default) (see [5])

    • ’qn1’ (see [6])

    • ’qn2’ (see [6])

    See here for more information.

  • likelihood (string, optional) –

    Likelihood approximation method. Should be one of (alphabetical ordering):

    • ’da’ (default) (see here)

    • ’Erlang’ (default) (see here)

    • ’expm’ (see here)

    • ’gwa’ (see here)

    • ’gwasa’ (see here)

    • ’ilt’ (see here)

    • ’oua’ (see here)

    • ’uniform’ (see here)

    The links point to the documentation of the relevant method in birdepy.probability(). The arguments associated with each of these methods may be used as a kwarg in birdepy.estimate() when likelihood is set to use the method.

  • laplace_method (string, optional) – Numerical inverse Laplace transform algorithm to use. Should be one of: ‘cme-talbot’ (default), ‘cme’, ‘euler’, ‘gaver’, ‘talbot’, ‘stehfest’, ‘dehoog’, ‘cme-mp’ or ‘gwr’.

  • lentz_eps (scalar, optional) – Termination threshold for Lentz algorithm computation of Laplace domain functions.

  • max_it (int, optional) – Maximum number of iterations of the algorithm to perform.

  • i_tol (scalar, optional) – Algorithm terminates when sum(abs(p(i) - p(i-1)) < i_tol where p(i) and p(i-1) are estimates corresponding to iteration i and i-1.

  • j_tol (scalar, optional) – States with expected number of upward transitions or expected number of downward transitions greater than j_tol are included in E steps.

  • h_tol (scalar, optional) – States with expected holding time greater than h_tol are included in E steps.

  • z_trunc (array_like, optional) – Truncation thresholds, i.e., minimum and maximum states of process considered. Array of real elements of size (2,) by default z_trunc=[z_min, z_max] where z_min=max(0, min(p_data) - 2*z_range) and z_max=max(p_data) + 2*z_range where z_range=max(p_data)-min(p_data).

Examples

Use birdepy.simulate.discrete() to simulate a sample path of the ‘Verhulst 2 (SIS)’ model:

import birdepy as bd
t_data = list(range(100))
p_data = bd.simulate.discrete([0.75, 0.25, 0.02, 1], 'Ricker', 10, t_data,
                              survival=True, seed=2021)

Estimate the parameter values from the simulated data using the ‘em’ framework with various technique and accelerator approaches:

for technique in ['expm', 'ilt', 'num']:
    for accelerator in ['cg', 'none', 'Lange', 'qn1', 'qn2']:
        tic = time.time()
        est = bd.estimate(t_data, p_data, [0.5, 0.5, 0.05], [[1e-6,1], [1e-6,1], [1e-6, 0.1]],
                          framework='em', technique=technique, accelerator=accelerator,
                          model='Ricker', idx_known_p=[3], known_p=[1])

The constraint con={'type': 'ineq', 'fun': lambda p: p[0]-p[1]} ensures that p[0] > p[1] (i.e., rate of spread greater than recovery rate).

A RuntimeWarning associated with SciPy’s minimize() function may appear, this can be ignored.

Notes

Estimation algorithms and models are also described in [1]. If you use this function for published work, then please cite this paper.

For a text book treatment on the theory of birth-and-death processes see [2].

See also

birdepy.estimate() birdepy.probability() birdepy.forecast()

birdepy.simulate.discrete() birdepy.simulate.continuous()

References

[View birdepy.estimate(framework=’em’) source code]

birdepy.estimate(framework=’lse’)

birdepy.interface_lse.discrete_est_lse(data, squares, model, b_rate, d_rate, z_trunc, idx_known_p, known_p, p0, p_bounds, con, opt_method, options)

Parameter estimation for discretely observed continuous-time birth-and-death processes using least squares estimation. See here for more information.

To use this function call birdepy.estimate() with framework set to ‘lse’:

bd.estimate(t_data, p_data, p0, p_bounds, framework='lse', squares='fm', z_trunc=())

The parameters associated with this framework (listed below) can be accessed using kwargs in birdepy.estimate(). See documentation of birdepy.estimate() for the main arguments.

Parameters:

  • squares (string, optional) –

    Squared error from the expected value approximation method. Should be one of (alphabetical ordering):

    • ’expm’

    • ’fm’ (default)

    • ’gwa’

  • z_trunc (array_like, optional) – Truncation thresholds, i.e., minimum and maximum states of process considered. Array of real elements of size (2,) by default z_trunc=[z_min, z_max] where z_min=max(0, min(p_data) - 2*z_range) and z_max=max(p_data) + 2*z_range where z_range=max(p_data)-min(p_data). Only applicable to squares method ‘expm’.

Examples

Simulate a sample path and estimate the parameters using the various likelihood approximation methods:

import birdepy as bd
t_data = list(range(100))
p_data = bd.simulate.discrete([0.75, 0.25, 0.02, 1], 'Ricker', 10, t_data,
                              survival=True, seed=2021)
for squares in ['expm', 'fm', 'gwa']:
    est = bd.estimate(t_data, p_data, [0.5, 0.5, 0.05], [[1e-6,1], [1e-6,1], [1e-6, 0.1]],
                      framework='lse', model='Ricker', idx_known_p=[3], known_p=[1],
                      squares=squares)
    print('lse estimate using', squares, 'is', est.p, 'computed in ', est.compute_time, 'seconds.')

Outputs:

lse estimate using expm is [0.7879591925854611, 0.26289368236374644, 0.02000871498805996] computed in  1.2579967975616455 seconds.
lse estimate using fm is [0.7941603523732229, 0.2766621569715867, 0.019363240909074483] computed in  5.483000755310059 seconds.
lse estimate using gwa is [0.7024952317382023, 0.20563650045779058, 0.022598851311981704] computed in  0.09800028800964355 seconds.

The constraint con={'type': 'ineq', 'fun': lambda p: p[0]-p[1]} ensures that p[0] > p[1] (i.e., rate of spread greater than recovery rate).

Notes

Estimation algorithms and models are also described in [1]. If you use this function for published work, then please cite this paper.

For a text book treatment on the theory of birth-and-death processes see [2].

See also

birdepy.estimate() birdepy.probability() birdepy.forecast()

birdepy.simulate.discrete() birdepy.simulate.continuous()

References

[View birdepy.estimate(framework=’lse’) source code]

birdepy.estimate(scheme=’continuous’)

birdepy.interface_dnm.continuous_est_dnm(t_data, p_data, p0, b_rate, d_rate, p_bounds, con, known_p, idx_known_p, opt_method, options)

Parameter estimation for continuously observed continuous-time birth-and-death processes using direct numerical maximization of the likelihood. See here for more information.

To use this function call birdepy.estimate() with scheme set to ‘continuous’:

bd.estimate(t_data, p_data, p0, p_bounds, scheme='continuous')

Examples

Simulate a continuous sample path and estimate the parameters.

Import BirDePy:

import birdepy as bd

Simulate some synthetic data:

t_data, p_data = bd.simulate.continuous([0.75, 0.25, 0.02, 1], 'Ricker', 10,
                                        100, survival=True, seed=2021)

Estimate:

est = bd.estimate(t_data, p_data, [0.5, 0.5, 0.05], [[0,1], [0,1], [0, 0.1]],
                  model='Ricker', idx_known_p=[3], known_p=[1], scheme='continuous')
print(est.p)

Outputs:

[0.7603171062895576, 0.2514810854871476, 0.020294342655751033]

Notes

Estimation algorithms and models are also described in [1]. If you use this function for published work, then please cite this paper.

For a text book treatment on the theory of birth-and-death processes see [2].

See also birdepy.estimate() birdepy.probability() birdepy.forecast()

birdepy.simulate.discrete() birdepy.simulate.continuous()

References

[View birdepy.estimate(scheme=’dnm’) source code]