How to make a parameter 2d in hierarchical model?

I want to model a hierarchical normal distribution with a 2d location parameter matrix mu and it’s regression parameters mu0 and mu1 which are hierarchical with hyperpriors alpha0 and alpha1. I also implement covariates COV to mu which have a specific length D.

See the code below:

data {
  int<lower=1> N;  \\sample size
  int<lower=1> J;   \\sample number
  int<lower=1> D;   \\3rd dimension
  vector[N] y[J,D];    \\sample 3d matrix
  vector[D] COV;    \\covariate for dimension

parameters {
  real<lower=0> sigma[J];
  real mu0[J];
  real mu1[J];

  real alpha0;
  real alpha1;

transformed parameters {
  real mu[J,D]
  mu = mu0 + mu1 * COV

model {
  mu0 ~ normal(alpha0,10)
  mu1 ~ normal(alpha1,10)
  for (dd in 1:D) {
    for (jj in 1:J) {
      y[jj,dd] ~ normal(mu[jj,dd],sigma)

How would I define mu in the transformed parameters block and then implement it in the model? Is this method correct?

I think you should have a loop over the dimensions of mu to cover all entries of the array, e.g. a nested loop:

for (j in 1:J) {
    for (d in 1:D) {
        mu[j,d] = your_function(mu0[j], mu1[j], COV[d]);

where your_function can be something like mu0[j] + mu1[j]*COV[d], but you have to check what produces the right entries for the matrix mu. Is that what you are looking for?