Vectors of vector parameters in Turing.jl?

General Usage
1

I want to run a Bayesian model in which there’s an overall US average behavior, which is the mean for a set of regional behaviors, and that is the mean for all the states in the region… the behavior is defined in terms of a basis expansion which has a certain number of coefficients… let’s say it’s 10 coefficients.

Defining the overall US average is easy:

nationcoefs ~ MvNormal(10,2.0)

Now, suppose I have defined 8 regions, and want to do something like:

for i in 1:8
   regioncoefs[i] ~ MvNormal(nationcoefs,scale)
end

And then for each state want something like:

for i in 1:51
   statecoefs[region[i]] ~ MvNormal(regioncoefs[region[i]],scale2)
end

and then I’ll set up a likelihood based on the expansions and the observed data.

The issue is how do I define/declare/inform Turing that “regioncoefs” is a certain sized vector of vectors, and similarly for “statecoefs”?

Can I just declare something like:

regioncoefs = [zeros(Float64,k) for i in 1:nregion]

before I do the ~ statements?

2

Ok, I’ve tried to do something like this, but am getting an error:

    regioncoefs = fill(nationcoefs,nregions)
    for i in 1:length(regioncoefs)
        regioncoefs[i] ~ MvNormal(nationcoefs,2.0)
    end

The error is:

ERROR: ArgumentError: Converting an instance of ReverseDiff.TrackedReal{Float64, Float64, Nothing} to Float64 is not defined. Please use `ReverseDiff.value` instead.

And I’m not sure it’s related to this code, but it seems to be that others have had similar issues with mutation based code.

I can’t figure out any other way to express the idea that there’s a vector of vectors and each vector in the collection should be distributed as some Distribution type. I’m not sure if something like

thing = [ foo ~ MvNormal(bar,scale) for i in 1:n]

is legal syntax / meaningful. but I’ll give that a try too

3

The error I was getting was because I was calling

sample(model,NUTS(0.8,500))

instead of

sample(model,NUTS(500,0.8))

So, I can basically say that yes it work to do the following:

    nationcoefs ~ MvNormal(zeros(Float64,nc),2.0)
    regioncoefs = fill(nationcoefs,nregions)
    for i in 1:length(regioncoefs)
        regioncoefs[i] ~ MvNormal(nationcoefs,2.0)
    end

It’s important to avoid doing something like fill(0.0,nregions) because it screws up the type system when Turing starts to use ReverseDiff types or ForwardDiff types etc.