Prototype binomial obs model

[seabbs]

@SamuelBrand1 any thoughts. The issue with the current approach is N needs to be a vector and the prior function has no way to support this. One option is a data unpacking generic method (empty) that can be specialised and other options are to doo with changing what is passed in in the first place

[github-actions]

Try this Pull Request!

Open Julia and type:

import Pkg
Pkg.activate(temp=true)
Pkg.add(url="https://github.com/CDCgov/Rt-without-renewal", rev="binomial-obs", subdir="EpiAware")
using EpiAware

[SamuelBrand1]

I think ultimately this goes back to #161 because we do want to lay out a long term data structure plan. Then the observed cases and the observed $N$ can both be equally handled here.

[SamuelBrand1]

@SamuelBrand1 any thoughts. The issue with the current approach is N needs to be a vector and the prior function has no way to support this. One option is a data unpacking generic method (empty) that can be specialised and other options are to doo with changing what is passed in in the first place

Total brainstorm (no wrong answer) first thought.

We can splat both positional and kwargs so a pattern that could work:

for i in eachindex(Y_t)
        y_t_defined[i + diff_t] ~ observation_error(obs_model, Y_t[i]...; priors...)
    end

This would require Y_t to be either a Vector{...} or a Vector{Vector{...}} in the vector of vector case each element would be [prob, N] which would get splatted out to the observation_error method dispatching on BinomialError. This does however seem to shift the problem one step up into defining Y_t. This pattern would work with BetaBinomial too I think by having the two positionals and the dispersion param.

Also padding would have to be done at the level down; which might make sense here because padding something on [0,1] is different

[SamuelBrand1]

In fact given our transform expectation utility maybe just get rid of padding (different issue)

[seabbs]

We can splat both positional and kwargs so a pattern that could work:

This is close to what we have and doesn't help because it still means N (which we are assuming is data and so not in Y_t which is model generated) is a vector?

[SamuelBrand1]

We can splat both positional and kwargs so a pattern that could work:

This is close to what we have and doesn't help because it still means N (which we are assuming is data and so not in Y_t which is model generated) is a vector?

N being a vector is a problem for the generate from prior approach unless we have some thing where we also pass in the index, e.g.

function observation_error(obs_model::BinomialError, Y_t, t ; N)
    return Binomial(N[t], Y_t)
end

I'm conceptually ok with N being generated by a model as a vector since we can condition on that vector. That ends up with us having the capability to generate N into the future (e.g. testing shifts etc).

[seabbs]

yes my point was your approach didn't seem to work?

[SamuelBrand1]

We can splat both positional and kwargs so a pattern that could work:

This is close to what we have and doesn't help because it still means N (which we are assuming is data and so not in Y_t which is model generated) is a vector?

The maths means we want/need $(N_t)_{t\geq0}$. That mean that the minimum data per time point is $(y_t, ~ N_t)$. That suggests that to support this we probably longer term want

...options are to doo with changing what is passed in in the first place

Because $N_t$ is as good a bit of data as $y_t$ so it seems a bit weird to have different ways of conditioning on one vs the other.

[SamuelBrand1]

More chucking ideas out there, instead of inelegant double splatting

for i in eachindex(Y_t)
        y_t_defined[i + diff_t] ~ observation_error(obs_model, Y_t[i], priors...)
 end

could be (yet another) @submodel

for i in eachindex(Y_t)
        @submodel generated_y_t = observation_error(obs_model, y_t_defined[i + diff_t] , Y_t[i], priors...)
 end

This gives more flexibility to pass in the data as a Vector{SomeDataObj}

e.g.

@model function observation_error(obs_model::BinomialError, y_t, Y_t)
	cases, N = y_t
	cases ~ Binomial(N, Y_t)
	return cases
end

I realise this is not cost free