Non-centered parameterization with boundaries for multiple coefficients

Hi, I would like to know if it is possible to declare bounded random effects with non-centered parameterization for multiple coefficients.

I read a very elegant solution here. This is great for one coefficient.

I am trying to do this for multiple coefficients. Currently my code goes like this:

data {
  int<lower=0> I; // number of groups
  int<lower=0> K; // number of coefficients
  int<lower=0> R; // number of  group-level covariates
  matrix[I, R] Z; // group-level covariates (with intercept)

  matrix[R, K] beta; //mean-slopes for all groups
  row_vector<lower=0>[K] sig; // scales for group-level coefficients
  matrix[I, K] eta;  // std_normal errors for group-level coefficients

transformed parameters{
  matrix[I, K] zeta; // group-level coefficients

  for (k in 1:K)
    zeta[,k] = Z*beta[,k] + eta[,k]*sig[k];

When the group-level coefficients are unbounded, I just need to specify the the number of variables and this code runs quickly and well.

I don’t know how to code it up so that it can handle the boundaries in my situation, even if the boundaries are the same for all the coefficients. I understand that if the number of coefficients (K) is fixed, I can just code up K different zeta’s brute-force. But I would like to know if there is a more flexible solution.

If one can offer any suggestion, I would much appreciate it.