Investigate rand calls in EpiAware

[SamuelBrand1]

One thing that became clear from trying to resolve #340 and doing #418 is that somewhere in an EpiAware model inference we are getting rand calls that triggered an error at large values of the mean.

However, its actually quite hard to recreate a MFE of this as per TuringLang/Turing.jl#2301

For example

using Turing, Distributions

@model function rw_model(rw₀, n)
    ϵ_t ~ filldist(Normal(), n)
    rw = rw₀ .+ cumsum(ϵ_t)
    return rw
end

@model function data_model(ys, rw)
    for i in eachindex(rw)
        ys[i] ~ Poisson(exp(rw[i]))
    end
    return ys
end

@model function pois_rw(ys, rw₀)
    n = length(ys)
    @submodel rw = rw_model(rw₀, n)
    @submodel gen_ys = data_model(ys, rw)
    return rw, gen_ys
end

generative_mdl = pois_rw(fill(missing,10), 3.)
ys = generative_mdl()[2] .|> Int
inference_mdl = pois_rw(ys, 48.) #deliberately bad choice for rw_0

chn = sample(inference_mdl, NUTS(), 1000)

works without needing SafePoisson. Clearly in this example the problematic rand call is being avoided despite being roughly analogous to an EpiAware model.

We should watch out for this, because it might indicate a suboptimal decision in design somewhere that we could fix more directly (as opposed to patching our own custom distribution types).