Vectorized Count Model

[Devin Caughey]

Dear all,

I am fitting a binomial count model, where the counts are stored in R as a 3D array (group by question by time). The model runs fine in non-vectorized form, with the array of counts declared as follows:
  int<lower=0> RRY[G, Q, T]; // count of correct answers in group, by year
  int<lower=0> NNY[G, Q, T]; // count of responses in group, by year
and the response model written as follows:
  RRY[g, q, t] ~ binomial(NNY[g, q, t], prob[q, t, g]); // group counts
To speed up computation, I would like to vectorize binomial(), but this requires that I declare RRY and NNY as arrays of vectors, as in:
  vector[T] RRY[G, Q];
This causes an error, however, because binomial() requires that the response be integer type. I haven't come across any way to declare an integer vector, which makes me think that they are not supported. Is this correct? If so, is there any alternative way to vectorize this problem?

Thanks!
~Devin

[Bob Carpenter]

On 3/26/13 3:25 PM, Devin Caughey wrote:
>
> I am fitting a binomial count model, where the counts are stored in R as a 3D array (group by question by time). The
> model runs fine in non-vectorized form, with the array of counts declared as follows:
>
> int<lower=0> RRY[G, Q, T]; // count of correct answers in group, by year
> int<lower=0> NNY[G, Q, T]; // count of responses in group, by year
>
> and the response model written as follows:
>
> RRY[g, q, t] ~ binomial(NNY[g, q, t], prob[q, t, g]); // group counts
>
> To speed up computation, I would like to vectorize binomial(), but this requires that I declare RRY and NNY as arrays of
> vectors, as in:
>
> vector[T] RRY[G, Q];

Not quite. Here's the signature for binomial (from
the index):

binomial_log
(ints n, ints N, reals theta):real

An "ints" type means it can be an int or an int *array*,
and "reals" means it can be a real, real array, or real vector.

You're right that there are no vectors of int types.

There is a simple and then an extreme form of vectorization
you could do for this statement:

RRY[g, q, t] ~ binomial(NNY[g, q, t], prob[q, t, g]); // group counts

If you reorganize prob to be indexed prob[g,q,t], then you can
vectorize on t:

RRY[g,q] ~ binomial(NNY[g,q], prob[g,q]);

Or you could reorganize and vectorize on the longest dimension.

At a more extreme level, you can create one large G * Q * T array
and fill it in with your own indexing and vectorize that way.
It'll be more efficient, despite all the copying, but less readable.

RRY_vec int[G * Q * T];
NNY_vec int[G * Q * T];
prob_vec real[G * Q * T];
int pos;
pos <- 1;
for (g in 1:G) {
for (q in 1:Q) {
for (t in 1:T) {
RRY_vec[pos] <- RRY[g,q,t];
NNY_vec[pos] <- NNY[g,q,t];
prob_vec[pos] <- prob[q,t,g];
pos <- pos + 1;
}
}
}
RRY_vec ~ binomail(NNY_vec,prob_vec);

- Bob

[Devin Caughey]

Thank you, Bob! This worked like a charm.

Devin

--
You received this message because you are subscribed to a topic in the Google Groups "stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/JyHaPT6EgY0/unsubscribe?hl=en.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.

=========================

Devin Caughey

Assistant Professor
MIT Department of Political Science

Office: E53-463, (617) 324-4085
cau...@mit.edu
www.devincaughey.com

[Devin Caughey]

Hi Bob,

I actually haven't been able to make the fully vectorized version work. I had to split up and rewrite the template code you suggested in order to make it work. I declared separately prob_vec as a transformed parameter and put the following block of code in the transformed data block:

  int RRY_vec[T * Q * G];
  int NNY_vec[T * Q * G];

  int pos;
  pos <- 1;
  for (g in 1:G) {
    for (q in 1:Q) {
      for (t in 1:T) {

        RRY_vec[pos] <- RRY[t, q, g];
        NNY_vec[pos] <- NNY[t, q, g];

        pos <- pos + 1;
      }
    }
  }

I put the last line of code in the model block:

  RRY_vec ~ binomial(NNY_vec, prob_vec);

But even when I comment out the last line of code, I get a memory allocation error and R crashes. Since I'm using a toy dataset, the length of the vectors is only 640, so finding the memory for them is not the problem. I must be declaring something incorrectly. The exact error I get when I start compiling the Stan model is printed below:

TRANSLATING MODEL 'stan.code' FROM Stan CODE TO C++ CODE NOW.

COMPILING THE C++ CODE FOR MODEL 'stan.code' NOW.

SAMPLING FOR MODEL 'stan.code' NOW (CHAIN 1).

R(11778,0x7fff7df4c180) malloc: *** error for object 0x1365b8a00: incorrect checksum for freed object - object was probably modified after being freed.

*** set a breakpoint in malloc_error_break to debug

Process R abort trap: 6 at Tue Mar 26 17:20:02 2013

~Devin

On Mar 26, 2013, at 3:47 PM, Bob Carpenter <ca...@alias-i.com> wrote:

[Bob Carpenter]

You should not be getting a malloc error and crash no matter
what you do, so this is definitely a bug in our code somewhere.

Which version of RStan are you using? I'd like to try to isolate
the problem as an issue in RStan vs. Stan itself.

Could you send along the whole model? Feel free to send
it just to me (ca...@alias-i.com) if you don't want to share
with the world. Data would help to, but I can do without that.

You don't want to have to put these local variables in the
transformed parameter block where they'll be saved. You can
make them locals in one of two ways. You can declare them at
the top of the model block, as in

model {

int RRY_vec[T * Q * G];
int NNY_vec[T * Q * G];
int pos;

...compute values of RRY and NNY...
RRY_vec ~ binomial(NNY_vec, prob_vec);
}

or they need to be enclosed in their own curly braces within
the model:

{

int RRY_vec[T * Q * G];
int NNY_vec[T * Q * G];
int pos;

... compute values of RRY and NNY ...
RRY_vec ~ binomial(NNY_vec, prob_vec);
}

The curly braces define a block which determines the scope and local variables;
local variables can only be declared at the top of a block.

- Bob

On 3/26/13 5:37 PM, Devin Caughey wrote:
> Hi Bob,
>
> I actually haven't been able to make the fully vectorized version work. I had to split up and rewrite the template code
> you suggested in order to make it work.

> I declared separately prob_vec as a transformed parameter and put the following
> block of code in the transformed data block:

To make this

>> <mailto:stan-users+...@googlegroups.com>.

>> For more options, visit https://groups.google.com/groups/opt_out.
>>
>>
>

> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Andrew Gelman]

Hey--we're helping a political scientist!  Cool!

You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Devin Caughey]

Ha! Cool for me you mean. I just wrote an email out to the department here at MIT to tell them they should all start using Stan. Thanks for all your amazing work on this!

Devin

[Ben Goodrich]

On Tuesday, March 26, 2013 5:37:47 PM UTC-4, Devin Caughey wrote:

I put the last line of code in the model block:

  RRY_vec ~ binomial(NNY_vec, prob_vec);

But even when I comment out the last line of code, I get a memory allocation error and R crashes. Since I'm using a toy dataset, the length of the vectors is only 640, so finding the memory for them is not the problem. I must be declaring something incorrectly. The exact error I get when I start compiling the Stan model is printed below:

TRANSLATING MODEL 'stan.code' FROM Stan CODE TO C++ CODE NOW.

COMPILING THE C++ CODE FOR MODEL 'stan.code' NOW.

SAMPLING FOR MODEL 'stan.code' NOW (CHAIN 1).

R(11778,0x7fff7df4c180) malloc: *** error for object 0x1365b8a00: incorrect checksum for freed object - object was probably modified after being freed.

*** set a breakpoint in malloc_error_break to debug

Process R abort trap: 6 at Tue Mar 26 17:20:02 2013

When you encounter a problem like this when Stan is called from R, do it after first calling

set_cppo("debug")

and that will often reveal the problem. See also

https://groups.google.com/d/msg/stan-users/4gv3fNCqSNk/FXHj6fn0SC0J

Ben

[Devin Caughey]

Just to update everyone, Bob's second email identified the problem correctly: I had declared the local variables in the wrong place. Everything worked once I put them at the top of the model block.

~Devin

To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.