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?

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

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.