Supporting vector (/array) parameters in the monty DSL

[edknock]

Currently we support fitting array parameters but do not have support for them in the DSL, so here are some examples of how we might want to implement things. Focused on vectors but should extend to arrays fairly easily.

A simple example
Intuitively I think we would want to index similarly to odin, without indexing on the LHS

m <- monty_dsl({
  beta[] ~ Gamma(1, 1)
})

though I imagine we'd need to declare dimensions

m <- monty_dsl({
  beta[] ~ Gamma(1, 1)
  dim(beta) <- c(2)
})

Multiline equations
I think we'd support multiline equations in a similar way (particularly for boundary conditions, see dependency section)

m <- monty_dsl({
  beta[1] ~ Gamma(1, 1)
  beta[2] ~ Gamma(1, 2)
  dim(beta) <- c(2)
})

Fixed data
Passing in fixed data as arrays would be handy!

fixed <- list(beta_shape = ...,
              beta_rate = ...,
              n_beta = ...)
m <- monty_dsl({
  beta[] ~ Gamma(beta_shape[i], beta_rate[i])
  dim(beta) <- c(n_beta)
}, fixed = fixed)

Would we need to include dimensions of beta_shape and beta_rate? Could we do something like dim(beta) <- dim(beta_rate) or even n_beta <- length(beta_rate)?

We could have the same beta_shape but varying beta_rate

fixed <- list(beta_shape = ...,
              beta_rate = ...,
              n_beta = ...)
m <- monty_dsl({
  beta[] ~ Gamma(beta_shape, beta_rate[i])
  dim(beta) <- c(n_beta)
}, fixed = fixed)

Dependency
It would be nice to allow dependency - perhaps this needs some restriction e.g. beta[j] can depend on beta[i] if and only if j > i?

fixed <- list(beta_rate1 = ...,
              n_beta = ...)
m <- monty_dsl({
  beta[1] ~ Gamma(1, beta_rate1)
  beta[2:n_beta] ~ Gamma(1, scale = beta[i-1])
  dim(beta) <- c(n_beta)
}, fixed = fixed)

Joint/multivariate distributions
Ultimately we might want to support multivariate distributions (e.g. multivariate normal or Dirichlet)

fixed <- list(mu = ...,
              sigma = ...,
              n_beta = ...)
m <- monty_dsl({
  beta ~ MultivariateNormal(mu, sigma)
  dim(beta) <- c(n_beta)
}, fixed = fixed)