Why don't vectorized RNGs return vectors when given vector input?

Developers
1

I just saw these nice clarifications on the vectorized RNG return types

from @mitzimorris, but now I’m wondering why the return types are what they are. It seems strange that this isn’t allowed:

// assume mu and sigma declared like this in parameters
vector[N] mu;
real<lower=0> sigma;

// in generated quantities
vector[N] y_rep = normal_rng(mu, sigma);

Instead, even though mu is a vector I have to declare y_rep as a real array:

real y_rep[N] = normal_rng(mu, sigma);

I would have thought it would make sense for the return type to match the input type, i.e. if inputting arrays return arrays, if inputting vectors return vectors. Am I missing an obvious (or not obvious) reason why that would have been problematic? If not, are people open to that possibility?

1 Like
2

I have corrected the section of the user’s where this issue arose and added the following:

The vectorized form of the normal random number generator is used with the original predictors x and the model parameters alpha, beta , and sigma. The replicated data variable y_rep is declared to be the same size as the original data y , but instead of a vector type, it is declared to be an array of reals to match the return type of the function normal_rng . Because the vector and real array types have the same dimensions and layout, they can be plotted against one another and otherwise compared during downstream processing.

Hope this isn’t too hand-wavy.

2 Likes
3

I think your edits to the user’s guide are great. My question is really why this

choice was made in the language. That is, why was the choice made to always return an array even if that doesn’t match the input type (the vector mu in my example above).

4

perhaps @rybern or @bbbales2 can answer - is this a language thing or a math_lib thing?

5

I think we picked arrays because that’s the default type for having N of something. When you do multivariates, you get an array of vectors. When you do ints, you get an array of ints. When you do reals, you get an array of reals. In that way, the vectorized rngs are consistent on what they return.

It’s not obvious how to automatically pick a return type if mu and sigma were different types (row_vector and vector, for instance, which I think should work).

And then the why-doesn’t-the-assignment work goes back to strict Stan types.

2 Likes
6

Thanks @bbbales2!

7

Thanks, @bbbales2. That’s exactly what I was going to write. I’ll add that we can’t change the return types now without breaking backward compatibility.

Nevertheless, it does seem like it’d be more reasonable to have

normal_rng(vector, real)

return a vector even if the mixed cases are ambiguous.

My intent when designing the Stan types in the first place was that vector and matrix types would be used for matrix arithmetic and linear algebra and naturally multivariate operations like multi-normal location and covariance matrix parameters, but arrays would be used everywhere else. The thin end of the wedge was overloading the distributions to be neutral among row vectors, vectors, and arrays. The elementwise operations were just the icing on the cake when we allowed pow(vector) to return a vector.

Eigen, by way of example, cleanly separates arrays and matrices. Arrays are for elementwise operations and matrices for linear algebra and matrix arithmetic. It requires a lot of interconversion, which is made more efficient because they deal in expression templates, not dense base matrices or arrays.

8

Not to get sidetracked here, but maybe I missed when

this will be (and is this signature missing the power)? I’m still writing exp(vector .* log(vector))

9

That’s not actually exposed yet, but should be available in 2.24. We’ll be introducing an elementwise pow using the .^ operator, so you can use things like:

vector = vector .^ real
vector = vector .^ vector
vector = real .^ vector

Which you can also call using the pow function:

vector = pow(vector, real)
vector = pow(vector, vector)
vector = pow(real, vector)
3 Likes
10

Oops. I meant exp(vector).

1 Like
11

What happens with

vector .^ row_vector

and other mixed operations?

12

The inputs for the binary vectorisation (unless real) on the c++ side are required to have the same dimension and container type. So vector .^ row_vector would not be legal, but real[] .^ int[] would be legal.

However, this makes me think that we need to allow int[] arrays as inputs for all vector types. For a function like bessel_first_kind, which has the signature:

real bessel_first_kind(int v, real x)

If that’s going to be vectorised, the signatures would need to be:

real[] bessel_first_kind(int[] v, real[] x)
vector bessel_first_kind(int[] v, vector x)
row_vector bessel_first_kind(int[] v, row_vector x)
matrix bessel_first_kind(int[,] v, matrix x)

At the moment, only the first vectorised signature (with real[]) is possible, but I feel like the others should also be available

1 Like
13

That’d be great. We also just want int[] arguments because if someone writes { 1, 2, 3 } the reuslt is an int[3] type.

In general, I want to move to full covariant typing. That means because int is assignable to real, int[] should also be assignable to real[]. (Function application starts with assignment to the argument variables, so those should work the same way.)