Choosing appropriate data types for later access in loops

[glae...@gmail.com]

Dear Stan Users,

this is my first post and I just started in Stan (switching from JAGS), so please be patient with me ... :-)

My trial-by-trial (hierarchical) model processes two pieces of empirical data currently stored in R as: 

(1) M - a 3-dimensional matrix [3,nTrials,nSubjects]

(2) FB - a 2-dimensional matrix [nSubjects,nTrials]

In the model block I have 2 for-loops over subjects and trials. In the outer subjects-loop I need to access the then 2-dimensional matrix M, and the vector FB for that subject:

After some deliberation on the Stan data types, I think it is most appropriate to declare M as an array of matrices and FB as an array of vectors, so that I can access them in the following way:

data {

  ...

  matrix[128,3] M[337]; // 128 trials, 337 subjects

  vector[128] FB[337];

  vector[128] ONES[337]; // an array of vectors filled with 1's (for relating model predictions to the "virtual data" = all successes)

  ...

}

model {

  

  ... // declarations and sampling statements for hierarchical priors

  

  for (sub in 1:nSubjects) {

    matrix[128,3] m;

    vector[128] fb;

    vector[128] ones;

    vector[128] prob;

    matrix[128,3] a;

    m <- M[s];

    fb <- FB[s];

    ones <- ONES[s];

    

    ... // sample individual parameters from hierarchical priors

    a[1,1] <- 1/3;

    a[1,2] <- 1/3;

    a[1,3] <- 1/3;

    for ( t in 1:nTrials) {

      prob[t] <- dot_product(m[t],a[t]^d[sub]) / ( sum(a[t]^d[sub]) + machine_precision() ); // m = data, a = internal variable, d = model parameter

      ones[t] ~ bernoulli(prob[t]);

      ... // update internal variables with feedback (fb) etc.

    }

  }

}

I have 5 questions:

1. Does the Stan data type (array of matrices/vectors) seem like an appropriate choice?

2. Would the declaration in the data block work, if a re-ordered by data in R to be M[nSubjects,nTrials,3] and FB[nSubjects,nTrials]?

3. Do I have to pay attention to the vector orientation (row,column) in the dot_product(m[t],a[t]^d[sub]) expression?

4. Is there are more elegant way to assign 1/3 to the first row of matrix a other than assigning each element as above)?

5. Is ones[t] ~ bernoulli(prob[t]) good programming in Stan (it worked in JAGS) or should I rather use increment_log_prob(bernoulli_log(1,prob[t]))?

Thanks a lot for your insights.

Jan

[Bob Carpenter]

On Aug 1, 2014, at 11:55 AM, glae...@gmail.com wrote:

> Dear Stan Users,
>
> this is my first post and I just started in Stan (switching from JAGS), so please be patient with me ... :-)

We're very patient.

> My trial-by-trial (hierarchical) model processes two pieces of empirical data currently stored in R as:
> (1) M - a 3-dimensional matrix [3,nTrials,nSubjects]
> (2) FB - a 2-dimensional matrix [nSubjects,nTrials]
>
> In the model block I have 2 for-loops over subjects and trials. In the outer subjects-loop I need to access the then 2-dimensional matrix M, and the vector FB for that subject:
>
> After some deliberation on the Stan data types, I think it is most appropriate to declare M as an array of matrices and FB as an array of vectors, so that I can access them in the following way:
>
> data {
> ...
> matrix[128,3] M[337]; // 128 trials, 337 subjects
> vector[128] FB[337];
> vector[128] ONES[337]; // an array of vectors filled with 1's (for relating model predictions to the "virtual data" = all successes)

You should declare just a single vector[128] in transformed data,
fill it with ones using rep_vector(0,128), and then re-use it.

> ...
> }
> model {
>
> ... // declarations and sampling statements for hierarchical priors
>
> for (sub in 1:nSubjects) {
> matrix[128,3] m;
> vector[128] fb;
> vector[128] ones;
>
> vector[128] prob;
> matrix[128,3] a;

Stan is not BUGS. You don't need to declare all these local variables
to use them. This just incurs overhead in copying.

> m <- M[s];
> fb <- FB[s];
> ones <- ONES[s];
>
> ... // sample individual parameters from hierarchical priors
>
> a[1,1] <- 1/3;
> a[1,2] <- 1/3;
> a[1,3] <- 1/3;

Why are you only setting a[1] and then using a[t]? This will be a bug
because a[t] will be undefined for t > 1.

If you want each to be 1/3, then again, just declare a transformed
data variable and reuse it.

Again, Stan is not R. Evaluating 1/3 yields 0. Stan uses C++-style
arithmetic. You need to use 1/3.0 or 1.0/3 or 1.0/3.0. See the
manual for a discussion. If you do it as above, the compiler will
give you a warning about trying to use division on integers.

>
> for ( t in 1:nTrials) {
> prob[t] <- dot_product(m[t],a[t]^d[sub]) / ( sum(a[t]^d[sub]) + machine_precision() ); // m = data, a = internal variable, d = model parameter

If you're going to use m[t], then just declare M
as a two-dimensional array of vectors and use
M[s,t] in this expression. There is no copying involved
in just accessing an array.

I have no idea what you're trying to do by adding machine_precision to this.

> ones[t] ~ bernoulli(prob[t]);

Again, you don't need to assign prob[t] to a local variable. You can just
drop the expression inside as an argument.

Again, Stan is not R. Stan variables are typed. And bernoulli
requires integer arguments, not real arguments. I think what you
really want to do is this:

1 ~ bernoulli(prob[t]);

But it's probably better if you just vectorize and use

1 ~ bernoulli(prob);

> ... // update internal variables with feedback (fb) etc.
> }
> }
> }
>
> I have 5 questions:
> 1. Does the Stan data type (array of matrices/vectors) seem like an appropriate choice?

See above.

> 2. Would the declaration in the data block work, if a re-ordered by data in R to be M[nSubjects,nTrials,3] and FB[nSubjects,nTrials]?

As long as everything evaluates to the right type.

> 3. Do I have to pay attention to the vector orientation (row,column) in the dot_product(m[t],a[t]^d[sub]) expression?

No, but you might be able to simplify and speed things up by using
matrix arithmetic rather than repeated vector arithmetic and then
using a vectorized form of Bernoulli.

> 4. Is there are more elegant way to assign 1/3 to the first row of matrix a other than assigning each element as above)?

Yes, see above for why 1/3 won't work and what will work.

> 5. Is ones[t] ~ bernoulli(prob[t]) good programming in Stan (it worked in JAGS) or should I rather use increment_log_prob(bernoulli_log(1,prob[t]))?

Neither. Use 1 ~ bernoulli(prob); as described above.

-Bob

[glae...@gmail.com]

Dear Bob,

this is all great advice in an amazing response time ! :-) Thank you so much. (see inline for responses)

>   vector[128] ONES[337]; // an array of vectors filled with 1's (for relating model predictions to the "virtual data" = all successes)

You should declare just a single vector[128] in transformed data,
fill it with ones using rep_vector(0,128), and then re-use it.  

Thanks. I was looking for exactly that.

 

Why are you only setting a[1] and then using a[t]?  This will be a bug
because a[t] will be undefined for t > 1.

I actually left out part of the model code to shorten the original post. I am updating a[t+1] as part of the trial loop, so an updated version is available on the next iteration.

 

If you want each to be 1/3, then again, just declare a transformed
data variable and reuse it.

Will do. Thanks.

 

Again, Stan is not R.  Evaluating 1/3 yields 0.  Stan uses C++-style
arithmetic.  You need to use 1/3.0 or 1.0/3 or 1.0/3.0.  See the
manual for a discussion.  If you do it as above, the compiler will
give you a warning about trying to use division on integers.

Thank you so much for pointing that out. This probably saved me half a day of debugging :-)

 

>
>     for ( t in 1:nTrials) {
>       prob[t] <- dot_product(m[t],a[t]^d[sub]) / ( sum(a[t]^d[sub]) + machine_precision() ); // m = data, a = internal variable, d = model parameter

If you're going to use m[t], then just declare M
as a two-dimensional array of vectors and use
M[s,t] in this expression.  There is no copying involved
in just accessing an array.

Thanks again. I guess I will still need some time to get used to the Stan data types (and their potential). It's quite different form JAGS/BUGS.

 

I have no idea what you're trying to do by adding machine_precision to this.  

to avoid a devision-by-zero problem, because a[t] could potentially hold all zeros. I needed this in JAGS, so I thought I use something similar in Stan ...

 

>       ones[t] ~ bernoulli(prob[t]);

Again, you don't need to assign prob[t] to a local variable.  You can just
drop the expression inside as an argument.

Again, Stan is not R.  Stan variables are typed.  And bernoulli
requires integer arguments, not real arguments.  I think what you
really want to do is this:

  1 ~ bernoulli(prob[t]);

But it's probably better if you just vectorize and use

   1 ~ bernoulli(prob);

YES. Thank you. I wasn't sure if that was possible.

 

> 2. Would the declaration in the data block work, if a re-ordered by data in R to be M[nSubjects,nTrials,3] and FB[nSubjects,nTrials]?

As long as everything evaluates to the right type.

And I would see problems as warnings/errors at compile time?

 

> 3. Do I have to pay attention to the vector orientation (row,column) in the dot_product(m[t],a[t]^d[sub]) expression?

No, but you might be able to simplify and speed things up by using
matrix arithmetic rather than repeated vector arithmetic and then
using a vectorized form of Bernoulli.

But I am updating a[t+1] as part of the trial loop and I thought that requires a trial-by-trial updating of vectors, rather than a single step matrix arithmetic. But maybe I am mistaken there as well. 

Anyways, I am attaching the full model code with your suggestions already implemented (I think). If you have the time, I'd really appreciate if you took another look.

Thanks so much.

Jan

[Bob Carpenter]

You definitely don't want to be dividing by 0. That'll produce
either infinity or NaN depending on whether the numerator is
finite or 0.

>
> > ones[t] ~ bernoulli(prob[t]);
>
> Again, you don't need to assign prob[t] to a local variable. You can just
> drop the expression inside as an argument.
>
> Again, Stan is not R. Stan variables are typed. And bernoulli
> requires integer arguments, not real arguments. I think what you
> really want to do is this:
>
> 1 ~ bernoulli(prob[t]);
>
> But it's probably better if you just vectorize and use
>
> 1 ~ bernoulli(prob);
>
> YES. Thank you. I wasn't sure if that was possible.

Other than in assignment statements, a <- b, expressions are
not sensitive to their form as long as they are the right type.
In assignments, the left-hand side can only be an indexed
variable.

> > 2. Would the declaration in the data block work, if a re-ordered by data in R to be M[nSubjects,nTrials,3] and FB[nSubjects,nTrials]?
>
> As long as everything evaluates to the right type.
>
> And I would see problems as warnings/errors at compile time?

Yes.

>
> > 3. Do I have to pay attention to the vector orientation (row,column) in the dot_product(m[t],a[t]^d[sub]) expression?
>
> No, but you might be able to simplify and speed things up by using
> matrix arithmetic rather than repeated vector arithmetic and then
> using a vectorized form of Bernoulli.
>
> But I am updating a[t+1] as part of the trial loop and I thought that requires a trial-by-trial updating of vectors, rather than a single step matrix arithmetic. But maybe I am mistaken there as well.

Looking at the rest of the model, I don't see any way you
could do that with matrix arithmetic.

> Anyways, I am attaching the full model code with your suggestions already implemented (I think). If you have the time, I'd really appreciate if you took another look.

It must not compile and run yet, because it's "else if", not "elseif".

You don't want to declare a parameter with
constraint to lie in (0,1), as in:

real<lower=0,upper=1> r_mu;

and then have a truncation that's narrower:

r_mu ~ beta(A,B) T[L,U];

What you want to do is this:

real<lower=L,upper=U> r_mu;

and then you only need

r_mu ~ beta(A,B);

The truncation doesn't add anything to the log prob if L and
U are constants, and the constraint on r_mu makes sure the
model prob is non-zero for every value meeting the constraint.
Otherwise, if you had 0 < r_mu < L, you have a valid r_mu
for the constraint, but not in the model.

But you don't even need the sampling statement for r_mu,
because A = 1 and B = 1, so r_mu gets a uniform prior, which
is the default if you don't add anything in the model. Stan
is smart enough to just not do anything with the r_mu sampling
statement though other than error checking, so it won't be
a noticeable overhead.

Also, Stan tends to be more stable than JAGS or BUGS, so
you usually don't need to do things like .001 lower bounds
and .999 upper bounds.

Also, you can assign A <- 1 no problem. Integers will be promoted
to reals where required. It's just in the 1/3 case that nothing
is causing them to be promoted, so it's integer arithmetic.

As to putting assignments in the transformed parameters or as local
variables in the model is just a matter of whether they get saved
as parameters in the output.

- Bob

[glae...@gmail.com]

Dear Bob, thanks again for taking the time to advise me on my first Stan model, especially your comments about truncation. That's really helpful also for the future. I'll follow your suggestions and will run the model next week :-)

Cheers,
Jan