Public Documentation

Module for defining utility functions.

struct DirectSample <: AbstractEpiSamplingMethod

Sample directly from a Turing model.

Fields

• n_samples::Union{Nothing, Int64}: Number of samples from a model. If an integer is provided, the model is sampled n_samples times using Turing.Prior() returning an MCMChains. Chain object. If nothing, the model is sampled once returning a NamedTuple object of the sampled random variables along with generated quantities

struct HalfNormal{T<:Real} <: Distributions.Distribution{Distributions.Univariate, Distributions.Continuous}

Create a half-normal prior distribution with the specified mean.

Arguments:

• μ: The mean of the half-normal distribution.

Returns:

• A HalfNormal distribution with the specified mean.

Examples:

using EpiAware, Distributions

hn = HalfNormal(1.0)
# output
EpiAware.EpiAwareUtils.HalfNormal{Float64}(μ=1.0)

filter out all the values that are less than 0

rand(hn)
# output
0.4508533245229199

cdf(hn, 2)
# output
0.8894596502772643

quantile(hn, 0.5)
# output
0.8453475393951495

logpdf(hn, 2)
# output
-3.1111166111445083

mean(hn)
# output
1.0

var(hn)
# output
0.5707963267948966

Fields

• μ::Real

struct SafeNegativeBinomial{T<:Real} <: Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}

Create a Negative binomial distribution with the specified mean that avoids InExactError when the mean is too large.

Parameterisation:

We are using a mean and cluster factorization of the negative binomial distribution such that the variance to mean relationship is:

\[\sigma^2 = \mu + \alpha^2 \mu^2\]

The reason for this parameterisation is that at sufficiently large mean values (i.e. r > 1 / p) p is approximately equal to the standard fluctuation of the distribution, e.g. if p = 0.05 we expect typical fluctuations of samples from the negative binomial to be about 5% of the mean when the mean is notably larger than 20. Otherwise, we expect approximately Poisson noise. In our opinion, this parameterisation is useful for specifying the distribution in a way that is easier to reason on priors for p.

Arguments:

• r: The number of successes, although this can be extended to a continous number.
• p: Success rate.

Returns:

• A SafeNegativeBinomial distribution with the specified mean.

Examples:

using EpiAware, Distributions

bigμ = exp(48.0) #Large value of μ
σ² = bigμ + 0.05 * bigμ^2 #Large variance

# We can calculate the success rate from the mean to variance relationship
p = bigμ / σ²
r = bigμ * p / (1 - p)
d = SafeNegativeBinomial(r, p)
# output
EpiAware.EpiAwareUtils.SafeNegativeBinomial{Float64}(r=20.0, p=2.85032816548187e-20)

cdf(d, 100)
# output
0.0

logpdf(d, 100)
# output
-850.1397180331871

mean(d)
# output
7.016735912097631e20

var(d)
# output
2.4617291430060293e40

Fields

• r::Real

• p::Real

struct SafePoisson{T<:Real} <: Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}

Create a Poisson distribution with the specified mean that avoids InExactError when the mean is too large.

Arguments:

• λ: The mean of the Poisson distribution.

Returns:

• A SafePoisson distribution with the specified mean.

Examples:

using EpiAware, Distributions

bigλ = exp(48.0) #Large value of λ
d = SafePoisson(bigλ)
# output
EpiAware.EpiAwareUtils.SafePoisson{Float64}(λ=7.016735912097631e20)

cdf(d, 2)
# output
0.0

logpdf(d, 100)
# output
-7.016735912097631e20

mean(d)
# output
7.016735912097631e20

var(d)
# output
7.016735912097631e20

Fields

• λ::Real

accumulate_scan(
    acc_step::AbstractAccumulationStep,
    initial_state,
    ϵ_t
) -> Any

Apply the `accumulate` function to the `AbstractAccumulationStep` object.
This is effectively a optimised version of a for loop that applies the
`AbstractAccumulationStep` object to the input data in a single pass.

# Arguments
- `acc_step::AbstractAccumulationStep: The accumulation step function.
- `initial_state`: The initial state of the accumulation.
- `ϵ_t::AbstractVector{<:Real}`: The input data.

# Returns
- `state::AbstractVector{<:Real}`: The accumulated state as returned by the
`get_state` function from the output of the `accumulate` function.

# Examples
```julia
using EpiAware
struct TestStep <: AbstractAccumulationStep
    a::Float64
end

function (step::TestStep)(state, ϵ)
    new_state = step.a * ϵ
    return new_state
end

acc_step = TestStep(0.5)
initial_state = zeros(3)

accumulate_scan(acc_step, initial_state, [1.0, 2.0, 3.0])

function get_state(acc_step::TestStep, initial_state, state)
    return state
end

accumulate_scan(acc_step, initial_state, [1.0, 2.0, 3.0])
```

censored_cdf(
    dist::Distributions.Distribution;
    Δd,
    D,
    upper
) -> Any

Create a discrete probability cumulative distribution function (CDF) from a given distribution, assuming a uniform distribution over primary event times with censoring intervals of width Δd for both primary and secondary events.

NB: censored_cdf returns the non-truncated CDF, i.e. the CDF without conditioning on the secondary event occuring either before or after some time.

Arguments

• dist: The distribution from which to create the PMF.
• Δd: The step size for discretizing the domain. Default is 1.0.
• D: The upper bound of the domain. Must be greater than Δd. Default D = nothing

indicates that the distribution should be truncated at its upperth percentile rounded to nearest multiple of Δd.

Returns

• A vector representing the CDF with 0.0 appended at the beginning.

Raises

• AssertionError if the minimum value of dist is negative.
• AssertionError if Δd is not positive.
• AssertionError if D is shorter than Δd.
• AssertionError if D is not a multiple of Δd.

Examples

using Distributions
using EpiAware.EpiAwareUtils

dist = Exponential(1.0)

censored_cdf(dist; D = 10) |>
    p -> round.(p, digits=3)

# output
11-element Vector{Float64}:
 0.0
 0.368
 0.767
 0.914
 0.969
 0.988
 0.996
 0.998
 0.999
 1.0
 1.0

censored_pmf(
    dist::Distributions.Distribution,
    ::Val{:single_censored};
    primary_approximation_point,
    Δd,
    D
)

Create a discrete probability mass function (PMF) from a given distribution, assuming that the primary event happens at primary_approximation_point * Δd within an intial censoring interval. Common single-censoring approximations are primary_approximation_point = 0 (left-hand approximation), primary_approximation_point = 1 (right-hand) and primary_approximation_point = 0.5 (midpoint).

Arguments

• dist: The distribution from which to create the PMF.
• ::Val{:single_censored}: A dummy argument to dispatch to this method. The purpose of the Val

type argument is that to use single-censored approximation is an active decision.

• primary_approximation_point: A approximation point for the primary time in its censoring interval.

Default is 0.5 for midpoint approximation.

• Δd: The step size for discretizing the domain. Default is 1.0.
• D: The upper bound of the domain. Must be greater than Δd.

Returns

• A vector representing the PMF.

Raises:

• AssertionError if the minimum value of dist is negative.
• AssertionError if Δd is not positive.
• AssertionError if D is not greater than Δd.

Examples

using Distributions
using EpiAware.EpiAwareUtils

dist = Exponential(1.0)

censored_pmf(dist, Val(:single_censored); D = 10) |>
    p -> round.(p, digits=3)

# output
10-element Vector{Float64}:
 0.393
 0.383
 0.141
 0.052
 0.019
 0.007
 0.003
 0.001
 0.0
 0.0

censored_pmf(
    dist::Distributions.Distribution;
    Δd,
    D,
    upper
) -> Any

Create a discrete probability mass function (PMF) from a given distribution, assuming a uniform distribution over primary event times with censoring intervals of width Δd for both primary and secondary events. The CDF for the time from the left edge of the interval containing the primary event to the secondary event is created by direct numerical integration (quadrature) of the convolution of the CDF of dist with the uniform density on [0,Δd), using the censored_cdf function. The discrete PMF for double censored delays is then found using simple differencing on the CDF.

NB: censored_pmf returns a right-truncated PMF, i.e. the PMF conditioned on the secondary event occurring before or on the final secondary censoring window.

Arguments

• dist: The distribution from which to create the PMF.
• Δd: The step size for discretizing the domain. Default is 1.0.
• D: The upper bound of the domain. Must be greater than Δd. Default D = nothing

indicates that the distribution should be truncated at its upperth percentile rounded to nearest multiple of Δd.

Returns

• A vector representing the PMF.

Raises

• AssertionError if the minimum value of dist is negative.
• AssertionError if Δd is not positive.
• AssertionError if D is shorter than Δd.
• AssertionError if D is not a multiple of Δd.

Examples

using Distributions
using EpiAware.EpiAwareUtils

dist = Exponential(1.0)

censored_pmf(dist; D = 10) |>
    p -> round.(p, digits=3)

# output
10-element Vector{Float64}:
 0.368
 0.4
 0.147
 0.054
 0.02
 0.007
 0.003
 0.001
 0.0
 0.0

get_param_array(chn::MCMCChains.Chains) -> Any

Extract a parameter array from a Chains object chn that matches the shape of number of sample and chain pairs in chn.

Arguments

• chn::Chains: The Chains object containing the MCMC samples.

Returns

• param_array: An array of parameter samples, where each element corresponds to a single

MCMC sample as a NamedTuple.

Example

Sampling from a simple model which has both scalar and vector quantity random variables across 4 chains.

using Turing, MCMCChains, EpiAware

@model function testmodel()
    y ~ Normal()
end
mdl = testmodel()
chn = sample(mdl, Prior(), MCMCSerial(), 2, 1, progress=false)

A = get_param_array(chn)

get_state(
    acc_step::AbstractAccumulationStep,
    initial_state,
    state
) -> Any

Processes the output of the `accumulate` function to return the final state.

# Arguments
- `acc_step::AbstractAccumulationStep`: The accumulation step function.
- `initial_state`: The initial state of the accumulation.
- `state`: The output of the `accumulate` function.

# Returns
- `state`: The combination of the initial state and the last element of
  each accumulated state.

prefix_submodel(
    model::AbstractModel,
    fn::Function,
    prefix::String,
    kwargs...
) -> Any

Generate a submodel with an optional prefix. A lightweight wrapper around the @submodel macro from DynamicPPL.jl.

Arguments

• model::AbstractModel: The model to be used.
• fn::Function: The Turing @model function to be applied to the model.
• prefix::String: The prefix to be used. If the prefix is an empty string, the submodel is created without a prefix.

Returns

• submodel: The returns from the submodel are passed through.

Examples

using EpiAware, DynamicPPL

submodel = prefix_submodel(FixedIntercept(0.1), generate_latent, string(1), 2)
submodel
# output
Model{typeof(prefix_submodel), (:model, :fn, :prefix, Symbol("#splat#kwargs")), (), (), Tuple{FixedIntercept{Float64}, typeof(generate_latent), String, Tuple{Int64}}, Tuple{}, DefaultContext}(EpiAware.EpiAwareUtils.prefix_submodel, (model = FixedIntercept{Float64}(0.1), fn = EpiAware.EpiAwareBase.generate_latent, prefix = "1", var"#splat#kwargs" = (2,)), NamedTuple(), DefaultContext())

We can now draw a sample from the submodel.

rand(submodel)

scan(f::AbstractModel, init, xs) -> Tuple{Any, Any}

Apply f to each element of xs and accumulate the results.

f must be a callable on a sub-type of AbstractModel.

Design note

scan is being restricted to AbstractModel sub-types to ensure: 1. That compiler specialization is activated 2. Also avoids potential compiler overhead from specialisation on f<: Function.

Arguments

• f: A callable/functor that takes two arguments, carry and x, and returns a new carry and a result y.
• init: The initial value for the carry variable.
• xs: An iterable collection of elements.

Returns

• ys: An array containing the results of applying f to each element of xs.
• carry: The final value of the carry variable after processing all elements of xs.

Examples

```jldoctest using EpiAware

struct Adder <: EpiAwareBase.AbstractModel end function (a::Adder)(carry, x) carry + x, carry + x end

scan(Adder(), 0, 1:5) #output ([1, 3, 6, 10, 15], 15)

spread_draws(chn::MCMCChains.Chains) -> DataFrames.DataFrame

spread_draws(chn::Chains)

Converts a Chains object into a DataFrame in tidybayes format.

Arguments

• chn::Chains: The Chains object to be converted.

Returns

• df::DataFrame: The converted DataFrame.

∫F(dist, t, Δd) -> Any

Calculate the CDF of the random variable X + U where X has cumulative distriubtion function F and U is a uniform random variable on [0, Δd).

This is used in solving for censored CDFs and PMFs using numerical quadrature.