Poisson distributions and int variables

Modeling
1

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);
}
1 Like
2

You can create a local block in the parameters, transformed parameters, or model block,
and inside the local block you can create as many int variable as you want.
not sure that fits in with the model you’ve described.

3

If I create an int variable inside a local block, can it be referred to from other blocks in the model?

4

In short: you cannot directly estimate an integer parameter in Stan, and marginalizing over discrete parameters with unbounded support (such as those with Poisson priors) can be difficult when it’s possible at all.

I think a better approach may be to change to an ordinal regression response; since you’re trying to approximate your action potential threshold with a continuous value, it would make sense to give it a continuous prior. There’s a section in the manual on these.

4 Likes
5

Thanks Christopher!
I’m specifically interested in testing the model with two Poisson distributions for the two experimental conditions but a consistent mapping between the number of action potentials to the response. So unless I’m not understanding your suggestion, I’m not sure if ordinal regression will help here. Do you know if JAGS support integer parameters?
Thanks again,
-Matan

6

Yes.

For just the reason Christopher Peterson metnions, this is the only case where I wish Stan had integer parameters was for missing count data.

7

It is unfortunate that one cannot efficiently generate missing count data. I have such a nice negative binomial model that I am loath to give up. So I’m going to try and generate missing data from a log normal where I match the first two moments of the lognormal to the negative binomial. Totally adhoc but I’ll see how well it works.

8

One trick that might be useful is that you can decompose the neg. binomial into Poisson with a gamma-distributed mean and treat the gamma as the unknown. This can help because a Poisson sometimes lends itself to analytic marginalization.

If that doesn’t help, gamma tends to be the continuous counterpart to neg. binomial so you might be better off using it instead of the lognormal.

If you need further assistance with this I would suggest you start a new topic and provide more details about your actual problem.

1 Like