I am implementing the QR reparametrization of section 9.2. However, my model calls for an array of matrices.

```
data {
int<lower=0> N;
int T; // number of days in campaign
int K; // number of predictors
matrix[N, K] x[T];
vector[N] y[T];
}
```

So for the transformed data block I want to set up T Q_ast matrices with the modification that each QR parameter is a array of dimension T.

```
for (t in 1:T){
Q_ast[t] = qr_Q(z[t])[, 1:K] * sqrt(N - 1);
R_ast[t] = qr_R(z[t])[1:K, ] / sqrt(N - 1);
R_ast_inverse[t] = inverse(R_ast[t]);
}
```

However, I am not sure how I would recover beta using the formula from the manual.

`beta = R_ast_inverse * theta`

since in my case both R_ast_inverse and theta are arrays of dimension T, I will end up with T sets of betas. Do I take the average?