Hey experts,

I’m trying to model behavioural data using a Poisson distribution.

Data from each subject is given as a vector of 6 numbers, representing the number of times the subject responded with a certain response out of 6 possible ones (for example in the data [10, 30, 20, 10, 0, 10] the subject chose the second response 30 times). I assume in my model that responses are chosen as a function of the number of ‘action potentials’ of a hypothetical neuron. An example participant might respond ‘1’ if this number is less than 3, respond ‘2’ if it is in [3,5), respond ‘3’ if it is in [5,7), and so on.

[edit: it is actually the case that for each subject I have two such vectors - one for each experimental condition. I assume that the criteria for responses is the same between the two conditions, and all that changes is the underlying Poisson distribution. I omitted this part from the original post to keep things short].

My approach was to fit a lambda variable per subject, parametrizing the Poisson distribution from which number of '‘action potentials’ are sampled, as well as 5 thresholds for the different responses (in the example above these would be [3, 5, 7, …]). I then compute the probabilities of falling inside each bin using poisson_cdf and model the data using a multinomial distributions with these probabilities and the total number of trials.

My problem is that *Stan won’t allow int variables, yet the function poisson_cdf expects an integer for the ‘N’ (in my case threshold) variable. Importantly, in my case this N is a variable that I want to approximate.*

Many thanks (model follows),

-Matan

```
data {
int<lower=0> nratings; // number of possible decision x confidence responses
int<lower=0> counts[nratings]; //
int<lower=0> ntrials;
}
parameters {
real<lower=0> firingRate;
vector[nratings-1] criterion_change;
}
transformed parameters {
vector[nratings-1] criterion;
vector[nratings] p;
criterion[1] = criterion_change[1];
for (i in 2:nratings-1)
criterion[i] = criterion[i-1]+criterion_change[i];
p[1] = poisson_cdf(criterion[1], firingRate);
for (i in 2:nratings)
p[i] = poisson_cdf(criterion[i], firingRate) - poisson_cdf(criterion[i-1], firingRate);
p[nratings] = 1 - poisson_cdf(criterion[nratings-1], firingRate);
}
model {
firingRate ~ gamma(1.0,1.0);
for (i in 1:nratings-1)
criterion_change[i] ~ gamma(1.0,1.0);
counts~ multinomial(p);
}
```