Simulate a Stochastic Birth-Death Process (Exponential Growth)

Simulates a stochastic birth-death process over multiple time steps, which models population dynamics in an environment where both birth and death rates can vary over time. This process follows an exponential growth model, where population size fluctuates due to random birth and death events.

Usage

sim_stochastic_exponential(n0, lambda, mu, steps, delta_t)

Arguments

n0: Numeric value indicating the initial population size at time 0. This must be a non-negative integer (>= 0).
lambda: Numeric vector or scalar specifying the birth rates at each time step. If a scalar is provided, it is assumed to be constant for all time steps. If a vector is provided, it should be of length steps, specifying the birth rate at each step.
mu: Numeric vector or scalar specifying the death rates at each time step. Like lambda, if a scalar is provided, it is assumed to be constant for all time steps. If a vector is provided, it should be of length steps, specifying the death rate at each step.
steps: Integer specifying the number of time steps to simulate. Must be a non-negative integer.
delta_t: Numeric vector or scalar specifying the time step sizes. If a scalar is provided, it is assumed to be constant for all steps. If a vector is provided, it should be of length steps, specifying the time step size at each step.

Value

A tibble with the following columns:

time: The cumulative time at each step, starting from 0.
count: The population size at each time step, calculated after applying the stochastic birth-death process.
group: A factor that groups contiguous intervals where the birth and death rates are constant (i.e., lambda and mu remain the same over consecutive time steps).

Details

This function simulates a population following an exponential growth model, where the population size changes due to random births and deaths. The rates of birth (lambda) and death (mu) may vary over time. At each time step, the function computes the expected number of births and deaths based on the rates lambda and mu. A random draw is then made to determine whether a birth or death event occurs. The population size is updated accordingly.

The simulation assumes that the population size is never negative, so if a death event causes the population to drop to zero, it will remain zero thereafter.

If steps is set to 0, the function returns the initial population size n0 immediately, without any simulation.

Examples

# Simulate a stochastic birth-death process with constant rates
sim_stochastic_exponential(n0 = 10, lambda = 0.2, mu = 0.1, steps = 10, delta_t = 1)
#> # A tibble: 11 × 3
#>     time count group
#>    <dbl> <dbl> <dbl>
#>  1     0    10     0
#>  2     1    16     0
#>  3     2    17     0
#>  4     3    17     0
#>  5     4    18     0
#>  6     5    22     0
#>  7     6    30     0
#>  8     7    31     0
#>  9     8    37     0
#> 10     9    35     0
#> 11    10    40     0

# Simulate with varying rates
lambda_rates <- c(rep(0.3, 5), rep(0.1, 5), rep(0.9, 5))
mu_rates <- c(rep(0.1, 5), rep(0.15, 5), rep(0.3, 5))
sim_stochastic_exponential(n0 = 10, lambda = lambda_rates, mu = mu_rates, steps = 15, delta_t = .1)
#> # A tibble: 16 × 3
#>     time count group
#>    <dbl> <dbl> <dbl>
#>  1   0      10     0
#>  2   0.1    11     0
#>  3   0.2    12     0
#>  4   0.3    12     0
#>  5   0.4    14     0
#>  6   0.5    15     0
#>  7   0.6    14     1
#>  8   0.7    15     1
#>  9   0.8    15     1
#> 10   0.9    16     1
#> 11   1      15     1
#> 12   1.1    15     2
#> 13   1.2    13     2
#> 14   1.3    15     2
#> 15   1.4    14     2
#> 16   1.5    15     2