edit: I’m wrong, see Bob’s stuff
The code isn’t actually multiplying two matrices:
y ~ normal(rows_dot_product(beta[jj] , x), sigma);
beta[jj] is a row_vector, and x is a matrix, it just dots the row beta[jj] by each row of x.
Should be the same answer as
beta[jj] * x' or
x * beta[jj]' (the normal here isn’t gonna care if you feed it a vector or a row_vector)