turing_model part 1: fixed-effects

[storopoli]

Ok, this is a big PR. 50% of the package was done in this PR.

Mainly it implements everything for specifying and sampling models using @formula macro that do not have random-effects, i.e. the (1 | group), (x1 | group) or (1 + x1 | group) terms inside the @formula.

Implemented likelihoods:

1. normal
2. student-t
3. bernoulli
4. poisson
5. negative binomial

Datasets in data/

I am also adding 3 datasets from stan-dev/rstanarm:

1. kidiq
2. wells
3. roaches

Their license is GPL-3. They are used extensively in tutorials(see storopoli/Bayesian-Julia and also Gelman & Hill (2007), Gelman et al. (2013) (the BDA), Gelman et al. (2020) (RoS).

I know that @rikhuijzer hates data being hard-coded into a package but they are very small and are used for tests and tutorials (not implement yet, but in the roadmap)

The turing_model docstring deserves extra attention.

Relates to #2.

@yebai feel free to review or ask for others to review.

[rikhuijzer]

I know that @rikhuijzer hates data being hard-coded into a package but they are very small and are used for tests and tutorials (not implement yet, but in the roadmap)

But the drawbacks from Artifacts are negligible. It could be something like

function dataset(name::AbstractString)
     path = ""
     if name == "roaches"
           path = joinpath(Artifacts"roaches", "roaches.csv")
    elseif ...
    end
    return read_data(path)
end

It's easier than it looks: https://pkgdocs.julialang.org/v1/artifacts/. It clutters the diff less, makes switching to a newer version of the dataset easier, you even have to handle paths slightly less yourself because you can just say joinpath(Artifact"roaches", "roaches.csv") and it makes adding larger datasets in the same way at a later moment possible.

[rikhuijzer]

Yes nice stuff 💯

[rikhuijzer]

Can this closure be moved into a separate function? It's very difficult to read now whether it's a function call or function definition. E.g.

function construct_normal_model(y, X, predictors, ...)
     return @model function normal_model(
           ...
     end
end

[phipsgabler]

Besides that, it might also be interesting to be able to access the model definition from outside for other purposes. I don't know how realistic that is, but e.g., to allow people to just look at and work with the interiors.

[storopoli]

Yes! I was thinking in a turing_code function to return or print the Turing underlying Turing code.

I think that we can also have that function to be the main one. turing_model would just call turing_code and eval the string as a macro. Is that possible?

[phipsgabler]

Well, ideally that would work in DPPL itself. We have been thinking about adding the original expression to the model struct and use it for printing, which would be pretty minor change but just hasn't been considered to be of high urgency.

However, this will then only work on the instantiated model object, not the result of the @model itself, which after all is just a bunch of method definitions of the model evaluator function.

I wouldn't do eval within the function. If necessary, you could have a top-level dict of names to code and evaluate those on the top level, but I don't really see the advantage.

[storopoli]

Ok, I will try to work on something during the week. This week is nasty because there are a lot of final committees and also tons of final semester assessments to grade. But it is the final week of the academic semester for me. I will add the turing_code suggestion in the comment down below /

[rikhuijzer]

Same as above. With more than 120 lines, this method is too long to be easily readable

[storopoli]

Things to do:

• Break up the huge turing_model into smaller functions. Maybe make family a struct instead of a string? Yes Normal() is the default. It is a Distribution type and Normal.
• Have turing_model call turing_code and eval it. turing_code should be exposed also.
• Testing with priors instead of my_prior
• Testing with f = @formula(...) instead of inside the turing_model
• Testing Chains stuff with only one chn variable

[storopoli]

Ready to review again. I could not implement easily turing_code function because I need to figure how to parse the Prior structs DefaultPrior CustomPrior to be displayed in the turing_code function. I think we should leave this for the 0.2.0 release or future releases.

The whole turing_model API for non-hierarchical models is done. I've created a custom type for the likelihood called Model. I had to be creative with the naming to avoid conflicts with the Distributions.jl types because we need them in the namespace for users to specify custom priors.