High-dimensional array subset to matrix?

[BillWalker]

Say I have 3-dimensional data, like:

my_data <- array(rnorm(2000),c(100,4,5))

I'm then iterating through this data in my model:

data {
  int<lower=1> I;
  int<lower=1> J;
  int<lower=1> K;


  real x[I,J,K];
  ...
}
parameters {
  real b[I,J];
  ...
}
model {
  ...
  for(i in 1:I) {
    y[i] ~ multi_normal(b[i] * x[i],sigma_y);
  }
  ...
}

In R, I would expect

b[i,] %*% x[i,,]

to give me the correct vector. In stan, if I understand correctly, this code won't work because x[i] is a two-dimensional array of reals, not a proper matrix. Can I coerce it to be a matrix? Alternately, could I make x a vector of matrixes?

Or should I just always flatten everything to two dimensions as most?

Thanks,

BW

[Bob Carpenter]

The simplest relevant declarations would be:

row_vector[J] b[I]; // b is an I-array of J-row-vectors

matrix[J,K] x[I]; // x is an I-array of (JxK)-matrices

though you could also define b with the same effect as:

matrix[I,J] b;

We followed MATLAB, not R/BUGS in our matrix operators,
so you can just use:

(b[i] * x[i])'

to produce a column vector (b[i] * x[i] is a row vector).

We are going to generalize multi_normal in an upcoming
release to allow row vectors (as well as 1D arrays)
of matching dimension as well as vectors as arguments.

You can also save a transpose by storing x in the opposite
order, with:

matrix[K,J] xt[I]; // xt[i] = x[i]'

vector[J] bt[I]; //

In this arrangement, (xt[i] * bt[i]) is a column vector
so you don't need the transpose. (You could even define
xt as transformed data so you could use the original data
coding.)

- Bob

[Andrew Gelman]

It might be a good idea in the manual to have an example like this, in which we format the same array in several different ways, so users can see the equivalencies.  We can then also explain why some parameterizations work more efficiently than others.

A

[Bob Carpenter]

On 9/11/12 4:56 PM, Andrew Gelman wrote:
> It might be a good idea in the manual to have an example like this, in which we format the same array in several
> different ways, so users can see the equivalencies.

There is already such a discussion. More than
one, in fact. Check out the section on vectorization
in the optimization chapter of the user's guide.

There should probably be a whole section in the user's
guide on vectors, row vectors, matrices and arrays.
One thing that's not adequately discussed with examples
is the return type for matrix operations. It's all
fully documented, but I think users coming from R who
are not used to reading function-signature doc may get
confused.

> We can then also explain why some parameterizations work more
> efficiently than others.

I can write up things like bernoulli vs. bernoulli_logit,
but I could use some help on the multivariate modeling side.

- Bob