How to return an array of vectors from a function?

Modeling
1

Here’s a short dummy example of what I want to do:

functions {
  vector[] sim(real r1, vector v1) {
    vector[2] v2 = to_vector({2.0, 4.0});
    // do stuff
    return({v1*r1, v2});
  }
}

transformed parameters {
  vector[2] v1 = to_vector({1.0, 2.0});
  vector[2] sim_result[2];
  
  sim_result = sim(2, v1);
  print(sim_result)
}

The function sim is supposed to return an array of vectors. It works fine when exposed to R with expose_stan_functions. But doesn’t compile when the function is called in the model or transformed parameters block. The error message is very long, but starts with “ERROR(s) during compilation: source code errors or compiler configuration errors!”.

Is it possible to return an array of vectors from a function or do i need to convert to something else first, like a matrix?

2

The error seems to be due to the {v1*r1, v2} part. It works when I code it as:

vector[] sim(real r1, vector v1) {
    vector[2] v2 = to_vector({2.0, 4.0});
    vector[2] res[2];
    // do stuff
    res[1, ] = v1*r1;
    res[2, ] = v2;
    return(res);
  }

I can do array of vectors assignments of the type array_of_vectors = {v1, v2, ...} in the model block. But that seems not to be allowed in the function block?

3

Assignment is the same everywhere. But the functions block may only contain function declarations and definitions.

Your function can be implemented mor simply as:

vector[] sim(real r1, vector v1) {
  return { r1 * v1, [2, 4]' };
}

Not sure why you want to mix these two things together, but that’s the most straightforward way to code it.

You don’t need parens around returns (Stan’s like C/C++/Java rather than R in this regard).

4

Thanks Bob. That syntax looked much nicer.

However, your example code doesn’t compile on my computer. Here a reproducible example:

functions {
  vector[] sim1(real r1, vector v1) {
    return { r1 * v1, [2, 4]' };
  }
  
  vector[] sim2(real r1, vector v1) {
    vector[2] res[2] = { r1 * v1, [2, 4]'};
    return res;
  }
  
  vector[] sim3(real r1, vector v1) {
    vector[2] res[2];
    res[1, ] = r1 * v1;
    res[2, ] = [2, 4]';
    return res;
  }
}

parameters {
  real mu;
}

model {
  vector[2] v2 = [2, 4]';
  vector[2] sim_result[2];
  
  sim_result = sim1(2, v2);

  print(sim_result);
  sim_result[1, ] ~ normal(mu, 2);
}

Calling sim1 or sim2 from the model block in the example code above result in compiling error on my computer. While sim3 works fine. This is a bit unexpected behaviour or maybe a bug.

5

Yeah, return { ... } is known to not work, but

vector[] sim1(real r1, vector v1) {
    vector[rows(v1)] out[2];
    out = { r1 * v1, [2, 4]' };
    return out;
  }

should work.

6

Is there an open issue for this? If so, could you include a link? I didn’t know this was a known issue.

If not, let’s create an issue.

7

I´m sorry, but that doesn´t work either. Only the sim3 version compiles on my computer

8

I’m not exactly sure what’s broken in that program in the current RStan 2.17.3 (presumably it’s also broken in the latest CmdStan as their stanc should be the same). But it works on the CmdStan develop branch, so presumably our next release will fix the issue.

So I couldn’t find a bug in how things work now with several variants of this.