3D Array multiplication and indexing

Modeling
1

I am working on converting my function in R into the Stan format and I am a little confused on how I should code my 3D array. I need to input data into my 3D array and it must be able to undergo matrix multiplication. Any feed back on how to reformat this function would be helpful, but particularly if I should be using real Nf[ , , ] or matrix[ , ] Nf[ ] for my 3D array.

functions {
do.function(int mesh, int time, int Q, vector dat, vector Fe, real cv_q, real cv, real sigma.p, real R0, real S0, real steep, matrix K) {
real Nf[mesh, time, Q];
real RR;
vector[time] Nrand;

Nf[,1,] = rep_matrix(dat, Q) ; 
  
 Nf[,1,] <- Nf[,1,] + normal_rng(rep_vector(0, time), cv_q * dat);
  if (Nf < 0) { Nf == 0 }; // make sure there are no negatives
  
  for (t in 2:time){
    for (q in 1:Q){
    real E;
    real Recruits;
    
     Nrand = normal_rng(0, sigma.p);
      
      Nf[, t, q] = K .* Nf[, t-1, q]  * dx;  // K is a mesh x mesh matrix (200x200)
      E = sum(Fe * Nf[, t-1, q]) * dx ;
    
    // BH equation
    Recruits = (4 * steep * exp(R0) * E) / ((S0 * (1 - steep)) + (E * (5 * steep - 1)));
      RR = exp(normal_rng(log(Recruits) - ((cv * log(Recruits))^2)/2, cv * log(Recruits) ));
      
      Nf[, t, q] = Nf[, t, q] + (RR+ Nrand) ; 
    } // end of Q loop

   if (Nf < 0) { Nf == 0 }; // make sure there are no negatives
} // end of T loop
} // end of fxn

And I am getting the following error @ the line with: Nf[,1,] = rep_matrix(dat, Q)

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Info: assignment operator <- deprecated in the Stan language; use = instead.
Dimension mismatch in assignment; variable name = Nf, type = real[ , ]; right-hand side type = matrix.
Illegal statement beginning with non-void expression parsed as
  Nf[:, 1, :]
2

What exactly needs to be used in matrix multiplication? Is it an arbitrary 2D slice of the 3D array along any of the three dimensions? Or is it the stack of 2D slices along a particular dimension. In the latter case create a 3D data structure as a 1D array of matrices, as you suggest with

The former case won’t be quite so straightforward.

3

In that case I believe I need to do the latter with the matrix format. I am multiplying a 200 x 200 matrix (K) by a slice of the 3D array(vector of 200) to get an output vector of 200 for that time point. The for loop uses the previous time point to update the current time point.

The 3D array is developed in the code by adding a matrix (rep_matrix(dat, Q)) into the 3D array as initial values before the for loop. Regardless of how I format the Nf 3D array, it is still giving the same error for indexing. Do you have any suggestions for initializing a 3D array with a matrix if it is in the matrix format ( matrix [ ,] Nf[ ])?

4

The matrix elements of this array will be indexed as N[1, , ]. It looks like maybe you are using N[ , 1, ] instead? You can think of this array as a one-dimensional array where each element is a matrix. So to pull out a matrix, you need to use the first indexing position.

6

After trying this method a couple of ways I believe my Nf 3D array must be declared as

real Nf[Q, mesh, time]

In order to input the data that I want. I have a vector of data that I want to use to initialize the first column (time) for each Q array.

vector[200] dat;

for ( j in 1:Q){
rep = normal_rng(rep_vector(0, mesh), cv_q * dat); // random no. generator w/ mean 0 and sigma = cv * data
  z = to_vector(rep);

Nf[j, ,1] = dat + z ; // initializes the column for each Q array
}

but I get the following error:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Dimension mismatch in assignment; variable name = Nf, type = real[ ]; right-hand side type = vector.
Illegal statement beginning with non-void expression parsed as
  Nf[j, :, 1]

How do I initialize my 3D array to have the same starting column (time =1) for each Q array?

7

You’ve declared Nf to be an array of reals, slices of which are also arrays. You can’t assign a vector to an array in Stan. You can assign a column vector to the columns of a matrix.

I don’t feel capable of giving more specific advice because I don’t really understand what constraints you’re operating under. For example:

I don’t understand what you mean by this or why it would be true.

I’m also unsure whether I understand what you mean by “Q array.”

Overall, my advice is to format and type your 3D data structure in accordance with what you need to do downstream. Do you need to treat slices of this array as matrices for matrix operations? Then use an array of matrices. Do you need to treat slices as vectors? Then use either an array of matrices (which will give you access to the row- and column-vectors) or as a 2D array of vectors. Once you have made a good choice based on what happens downstream, then figure out how to coerce your input to the right format to play nicely.

8

Thanks @jsocolar. This last comment helped me rethink how to problem solve.