Fwd: Stan sampling speed with large dataset

[Ben Goodrich]

On Tuesday, January 8, 2013 4:43:29 PM UTC-5, Ben Goodrich wrote:

On Tuesday, January 8, 2013 4:07:17 PM UTC-5, Ben Goodrich wrote:

One thing you can do as of v1.1.0 is replace

  for (s in 1:S) {
    row(e_beta,s) ~ normal(0,1);
  }

with just

e_beta ~ normal(0,1);

Actually, this doesn't parse in Stan v1.1.0. Never mind.

Am I correct that something like

parameters {
  matrix[3,2] e_beta;
}
model {
  e_beta ~ normal(0,1);  
}

is supposed to parse in v1.1.0? Instead, it throws

DIAGNOSTIC(S) FROM PARSER:
no matches for function name="normal_log"
    arg 0 type=matrix
    arg 1 type=int
    arg 2 type=int
unknown distribution=normal
Parser expecting: "}"

Ben

[Bob Carpenter]

Vectorization for distributions only works for:

scalars
1D arrays
vectors
row vectors

I'm going to move the discussion of this more up front,
but I don't really know a good place to put it. Right
now, it's in section 1 of the Continuous Probabilities chapter
(currently that's section 26.1).

At the very least, we should fix the error message
so that it no longer says:

> unknown distribution=normal

It's right in that there's no normal for those arguments, but
the distribution is known, so it's confusing.

- Bob

> --
>
>

[Ben Goodrich]

On Tuesday, January 8, 2013 5:14:38 PM UTC-5, Bob Carpenter wrote:

Vectorization for distributions only works for:

    scalars
    1D arrays
    vectors
    row vectors

But the C++ with the VectorView things would handle a matrix or an array of vectors, right? Do we just need to update the parser?

Ben

[Bob Carpenter]

Yes, if we want to allow matrices, that should do it.
But...

The reason we didn't is that we never defined what
arguments should be legal.

One simple approach is to provide a linearization for
each data type (for 1D arrays it's natural, for matrices,
column major, etc.) The problem is that it would allow
users to write things like this:

matrix[3,2] y;
matrix[2,3] mu;
real sigma;

y ~ normal(mu,sigma);

where we would wind up with:

y[1,1] ~ normal(mu[1,1],sigma);
y[2,1] ~ normal(mu[2,1],sigma);
y[3,1] ~ normal(mu[1,2],sigma);
...
y[3,2] ~ normal(mu[2,3],sigma);

This would also allow:

matrix[3,2] y;
vector[6] mu;
real y;

y ~ normal(mu,sigma);

We'd also talked about what would happen if we did something
like:

vector[3] y;
matrix[4,3] m;

m ~ normal(y,sigma);

We'd want to repeat the column y 4 times?

I think everyone decided all of this is just too confusing
and that if there were vector, matrix, etc. args, they
all had to have the same shape (allowing for 1D arrays,
row vectors and column vectors as the same shape, so you
could do this:

vector[3] y;
row_vector[3] mu;
real<lower=0> sigma[3];

y ~ normal(mu,sigma);

This is already allowed.

- Bob


On 1/8/13 5:26 PM, Ben Goodrich wrote:

> --
>
>

[Ben Goodrich]

On Tuesday, January 8, 2013 5:39:27 PM UTC-5, Bob Carpenter wrote:

I think everyone decided all of this is just too confusing
and that if there were vector, matrix, etc. args, they
all had to have the same shape (allowing for 1D arrays,
row vectors and column vectors as the same shape, so you
could do this:

    vector[3] y;
    row_vector[3] mu;
    real<lower=0> sigma[3];

    y ~ normal(mu,sigma);

This is already allowed.

I think that is fine to have the same shape rule, with the proviso that a scalar is always legal.

Ben

[Bob Carpenter]

I think that makes the most sense. I'll have to check on
where the error-checking code is.

Mismatched linearization is simplest to code and doc.

I think requiring matching dimensions on multi-dim
inputs makes most sense (and yes, always allow simple
scalars everywhere).

It would also make sense to me stack row vectors
(or line up column vectors) to create a matrix, but
nobody seemed to think that would be useful.

- Bob

> --
>
>