Hi Howard,

I think I have found a better solutions! We can make use of the `joint`

function which allows users to combine two or more different distributions into a single likelihood.

For example: if we assume that vector `a`

is length 10 and each element is distributed \sim N(\mu_1, 1) and vector `b`

is length 10 and each element is distributed \sim N(\mu_2, 1) we could have a model like:

```
mu_1 <- normal(0, 100)
mu_2 <- normal(0, 100)
x <- cbind(rnorm(10, 10, 1), rnorm(10, 2, 1))
distribution(x) <- joint(normal(mu_1, 1, dim = 10), normal(mu_2, 1, dim = 10))
```

Notes:

- The use of
`cbind`

ensures that the data has the same dimensions as the collections of variables produced by `joint`

.
- When we extend the distributions with
`dim`

all the variables are independent

After reading the documentation for `joint`

I am confident that this will definitely give you the likelihood you want.

All the best!