stan questions

[Andrew Gelman]

I'm using Stan for research!

1. For some reason the page numbers are sporadic, for example no page numbers between 60 and 70. (I actually noticed this because I was actually looking for a particular page, having used the table of contents.)

2. It would help to have a single pdf document with the Jags examples with Jags and Stan code. I could then use this as a template. As is, I'm doing a lot of guessing, trying to figure out what will work and what won't, with syntax.

3. Can we put multiple data or parameter declarations in a single command line? I tried the following but it did not work:

int[0,] J.1, J.2, J.3; // number of persons, between-plot treatments, within-plot treatments

I switched to the following which was fine:

int[0,] J.1; int[0,] J.2; int[0,] J.3; // number of persons, between-plot treatments, within-plot treatments

But the first form above would obv be convenient if it is possible.

4. I have a vector variable, real theta.1[J.1], and I give it a prior all at once:

theta.1 ~ normal (0, s.1);

That seems to work. But I also have a two-way array, real theta.1.3[J.1,J.3], and when I try to give it a prior all at once, it fails:

theta.1.3 ~ normal (0, s.1.3);

Could it work to allow such a call? It would be convenient for hierarchical models.

5. I had the following line:

for (i in 1:n)
y[i] ~ normal (mu + theta.1[i.1[i]] + theta.2[i.2[i]] + theta.3[i.3[i]] + theta.1.3[i.1[i],i.3[i]] + theta.2.3[i.2[i],i.3[i]], s.error);

Would it be possible to write this R-style with the looping implicit, like this?

y ~ normal (mu + theta.1[i.1] + theta.2[i.2] + theta.3[i.3] + theta.1.3[i.1,i.3] + theta.2.3[i.2,i.3], s.error);

If not, that's fine, I'm not saying that Stan needs every possible feature. I'm just curious.

Right now I'm still getting compilation errors, am working on debugging my model...

A

P.S. Stan's parser errors are awesome.

[Bob Carpenter]

On 8/10/12 6:27 AM, Andrew Gelman wrote:
> I'm using Stan for research!
>
> 1. For some reason the page numbers are sporadic, for example no page numbers between 60 and 70. (I actually noticed this because I was actually looking for a particular page, having used the table of contents.)
>
> 2. It would help to have a single pdf document with the Jags examples with Jags and Stan code. I could then use this as a template. As is, I'm doing a lot of guessing, trying to figure out what will work and what won't, with syntax.

You should at least glance through the manual.

But I understand that most users won't, so this is a good
idea. I'd rather NOT include the JAGS examples, as I think
it'll be confusing for people scanning.

I would have written a simple BUGS/JAGS to Stan conversion
guide, but it's rather complicated.

> 3. Can we put multiple data or parameter declarations in a single command line? I tried the following but it did not work:
>
> int[0,] J.1, J.2, J.3; // number of persons, between-plot treatments, within-plot treatments

It's on my to-do list, so you'll get lists like this
in the future. That, and I want to be able to allow

int J <- 10;

as well in blocks that allow variable definitions (all but
data and parameters).

You can't have "J.1" as a variable name -- no periods, only
underscores. As usual, I followed the usual (i.e., C++, Java,
etc. approach to variable names), not BUGS/R.

> I switched to the following which was fine:
>
> int[0,] J.1; int[0,] J.2; int[0,] J.3; // number of persons, between-plot treatments, within-plot treatments
>
> But the first form above would obv be convenient if it is possible.
>
> 4. I have a vector variable, real theta.1[J.1], and I give it a prior all at once:
>
> theta.1 ~ normal (0, s.1);
>
> That seems to work. But I also have a two-way array, real theta.1.3[J.1,J.3], and when I try to give it a prior all at once, it fails:
>
> theta.1.3 ~ normal (0, s.1.3);
>
> Could it work to allow such a call? It would be convenient for hierarchical models.

It's possible, but a bit of pain. We have to generalize
the vectorization, which in this case, requires a real
calculation to deal with arbitrary arrays.

Then we need to decide how it'll work. Is it OK to
allow

theta ~ normal(mu,sigma)

if theta is a 4x3-array, mu is a 2x6-vector and
sigma is a 12-vector? What about mixing a 4x3 array
and a 4x3 matrix?

If the shapes have to match, it's going to be significantly
harder to allow arrays.

Error checking would need a lot of work, as would the
vector wrappers.

> 5. I had the following line:
>
> for (i in 1:n)
> y[i] ~ normal (mu + theta.1[i.1[i]] + theta.2[i.2[i]] + theta.3[i.3[i]] + theta.1.3[i.1[i],i.3[i]] + theta.2.3[i.2[i],i.3[i]], s.error);
>
> Would it be possible to write this R-style with the looping implicit, like this?
>
> y ~ normal (mu + theta.1[i.1] + theta.2[i.2] + theta.3[i.3] + theta.1.3[i.1,i.3] + theta.2.3[i.2,i.3], s.error);
>
> If not, that's fine, I'm not saying that Stan needs every possible feature. I'm just curious.

I had to fire up R to see what its behavior is here.

> a = c(1,2,3)
> a[c(2,1,1,1)]
[1] 2 1 1 1

If v is a vector and s is a scalar, we already have
(v + s) doing the right thing in treating s like
a vector with duplication.

So the feature we need is just to have an array, vector or
matrix with an index consisting of an int array to return
a vector.

This is the easiest to implement of the requests :-)

>
> Right now I'm still getting compilation errors, am working on debugging my model...
>
> A
>
> P.S. Stan's parser errors are awesome.

There aren't any tests on the messages themselves, so let
me know (with model source) if you find an error message
you can't understand.

- Bob

[Andrew Gelman]

On Aug 10, 2012, at 6:23 PM, Bob Carpenter wrote:

>>
>> 4. I have a vector variable, real theta.1[J.1], and I give it a prior all at once:
>>
>> theta.1 ~ normal (0, s.1);
>>
>> That seems to work. But I also have a two-way array, real theta.1.3[J.1,J.3], and when I try to give it a prior all at once, it fails:
>>
>> theta.1.3 ~ normal (0, s.1.3);
>>
>> Could it work to allow such a call? It would be convenient for hierarchical models.
>
> It's possible, but a bit of pain. We have to generalize
> the vectorization, which in this case, requires a real
> calculation to deal with arbitrary arrays.
>
> Then we need to decide how it'll work. Is it OK to
> allow
>
> theta ~ normal(mu,sigma)
>
> if theta is a 4x3-array, mu is a 2x6-vector and
> sigma is a 12-vector? What about mixing a 4x3 array
> and a 4x3 matrix?

That would be fine with me. I was actualy assuming that mu and sigma would be scalars.

[Bob Carpenter]

mu is a vector when there's an underlying
regression.

Thinking about using matrices by flattening
them to vectors (and thus allowing an easy
extension of the VectorView we have), we
immediately run into order issues.

Eigen will flatten a 2x3 matrix using column-major
order:

(1,1), (2,1), (1,2), (2,2), (1,3), (2,3)

and a 2 x 3 matrix in C++'s row-major order would be

(1,1), (1,2), (1,3), (2,1), (2,2), (2,3)

This will confuse people if we do it and make the
code slow for one of Eigen or std::vector if we don't.

Let me think about this some more.

- Bob

[Andrew Gelman]

Yes, having mu as a vector would be delightful. I did not think that was possible, but even scalar mu, sigma would be cool because it would allow a quick specification of the distribution of error terms with constant variance. So if it's too much of a pain in the ass to implement anything too fancy, we could just allow scalars, or we could do some easy to implement compromise such as allowing either scalars or objects that are identical in shape to the left-hand-side object, but nothing in between. (That's something like what R does, I think.)

[Bob Carpenter]

On 8/10/12 5:29 PM, Andrew Gelman wrote:
> Yes, having mu as a vector would be delightful. I did not think that was possible,

Please, please, RTFM!

It's what I'm going to reply to the list to
these kinds of questions.

- Bob