Re: [stan-users] Can a Multivariate Poisson Distribution be implemented in Stan?

[Bob Carpenter]

You're right that you can't code the direct model in Stan because
it involves discrete parameters. If the counts aren't too high, it'd
be possible in theory to marginalize out the possible x values by
enumerating them.

Or you could build a more direct model where you model the counts
with some correlation. That is, you let lambda[1,...,K] be drawn
from a multivariate distribution directly. I am not familiar with
any of these models, so don't really know what the goal is.

We often use overdispersed models where you have a bit more latitutde
than in the Poisson (which has only a single parameter, which is both
mean and variance).

I'm a little confused by your model because you seem to have swapped
x and y from the description, with X in the Stan program being y in
the description. But otherwise it looks OK.

You can code the model a lot more efficiently.

This

> for (k in 1:K) {
> lambda[k] ~ gamma(0.1, 0.1);
> }
>

can just be

lambda ~ gamma(0.1, 0.1);

though I'd recommend more informative priors than this.

and

> for (n in 1:N) {
> for (m in 1:M) {
> lmbd <- 0;
> for (k in 1:K) lmbd <- lmbd + A[m, k] * lambda[k];
> X[n, m] ~ poisson(lmbd);
> }
> }

should be (with A being coded as a matrix and lambda as a vector,
as

for (m in 1:M) {
lmbd <- 0;
for (k in 1:K) if (A[m, k]) lmbd <- lmbd + lambda[k];
X[ , m] ~ poisson(lmbd);
}

Even better if you can transpose X and just use X[m] ~ ...
Less copying that way. This isn't the most efficient way to do this,
but it should be a whole lot faster than what you have above.
And if the X[, m] isn't working in this context (2.10 fixes some
of the edge cases), then X[1:N, m] will.

And if you want to move this all to a GLM setting, you can use poisson_log,
which takes the log of the mean as a parameter.

- Bob

> On Apr 19, 2016, at 5:55 PM, Aleksandar Dimitriev <dimitriev....@gmail.com> wrote:
>
>
> Let's say that I have a random data set of N observations of a k-dimensional vector that is drawn from a MVP distribution. I will try to briefly describe the approach, but it's probably better outlined here around slide 30: http://www.stat-athens.aueb.gr/~karlis/multivariate%20Poisson%20models.pdf
> So, the model is the following:
> Let Yn = (y1, ..., ym)' be one observation. The assumption is that Y = AX where X=(x1,...xk)' is a vector of univariate uncorrelated Poissons, and A is an indicator matrix such that each element of Y is a sum of particular elements of X. For example, for a trivariate case we may have:
> Y = (y1, y2, y3)', and if we want to model covariances between each y we can create a vector X = (x1, x2, x3, x12, x23, x13). Thus, the latter 3 terms would model the (positive) covariance. My goal is given a set of draws from this distribution, to infer the rates of each of the independent Poissons in X.
> The problem is that I can't seem to implement this in Stan because the Y are not observed and therefore should probably be parameters, but can't be since they are discrete, so the following code does not compile, of course.
> One thing that managed to work in Stan is another model: since each Y is a sum of certain X Poissons it is therefore distributed Poisson with a rate that is the sum of the rates of its Poissons, but I'm pretty sure this model is not equivalent to the previous one. This model in Stan is:
>
> data {
> int N; // number of observations
> int M; // M-variate Poisson
> int K; // M + M(M - 1) / 2, i.e. number of mean and covariance latent Poissons
> int A[M, K]; // indicator matrix
> int X[N, M]; // observed data
> }
>
> parameters {
> real<lower=0> lambda[K]; // rates of the latent Poissons
> }
>
> model {
> real lmbd;
> for (k in 1:K) {
> lambda[k] ~ gamma(0.1, 0.1);
> }
>
> for (n in 1:N) {
> for (m in 1:M) {
> lmbd <- 0;
> for (k in 1:K) lmbd <- lmbd + A[m, k] * lambda[k];
> X[n, m] ~ poisson(lmbd);
> }
> }
> }
>
> Also, the results from 2000 samples of this model had very low effective sizes and inference took a long time (~15min for a dataset of 1000 points of a trivariate Poisson). For reference, I generated a toy data with mean rates (lambda 1,2,3) of 10, 40 and 2, and covariance rates (lambda 12 13 and 23, or 4 5 and 6 in Stan output) of 50, 5 and 8, and I am a bit puzzled with the close-to-zero values of 1 3 and 5.
>
>
> Inference for Stan model: karlis.
> 4 chains, each with iter=1000; warmup=500; thin=1;
> post-warmup draws per chain=500, total post-warmup draws=2000.
> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
> lambda[1] 0.31 0.10 0.88 0.00 0.00 0.00 0.13 2.98 71 1.06
> lambda[2] 19.83 1.47 4.41 17.54 18.19 18.57 19.15 36.13 9 1.16
> lambda[3] 0.34 0.30 1.47 0.00 0.00 0.00 0.05 3.45 24 1.08
> lambda[4] 64.14 0.63 2.12 56.86 64.40 64.76 64.98 65.34 11 1.12
> lambda[5] 0.50 0.42 1.93 0.00 0.00 0.00 0.08 6.87 21 1.10
> lambda[6] 14.26 0.85 2.58 3.87 14.71 14.99 15.11 15.29 9 1.15
> lp__ 584414.88 0.20 2.02 584409.92 584413.77 584415.24 584416.32 584417.82 102 1.04
>
> Lastly, I saw that the bivariate case has been handled in Stan (https://groups.google.com/forum/#!newtopic/stan-users/stan-users/9mMsp1oB69g), but this is just hand-coding the bivariate PMF.
> I've also attached a traceplot of the 4 chains, does it indicate some problem with identifiability or is the likelihood just multi-modal?
>
> --
> 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <traceplot.png>

[andre.p...@googlemail.com]

Aleksandar,

I get a little confused of your Ys,Xs too Thus further I'm refering to pg. 33 of the pdf. 

I agree the observation of Y are discrete, but why X1, ..., X3 should be discrete,

when Y1,... are not? 

Isn't the matrix not just the sum of Poisson(s)? 

Andre

[Herra Huu]

As Aleksandar suspected, the current model doesn't correspond to MV Poisson distribution. In the current form it's just the product of marginal distributions, ie. defining joint density of form P(X1,X2,X3) = P(X1)P(X2)P(X3) and therefore missing the dependency structure.

The correct joint density (for M=3) is given in page 36. Basically what you have to do is:

take a sum over product of p(y1)p(y2)p(y3)p(y12)p(y13)p(y23), where each term is Poisson density, over all possible combinations of (y1,y2,y3,y12,y13,y23) such that:
x1 == y1+y12+y13
x2 == y2+y12+y23
x3 == y3+y13+y23

From there, it's easy to generalize for higher dimensions/different dependency structures, but it's also awfully computationally expensive if M/counts are high. Here is a simple (but not very efficient) example code for the M=3 case:
data generation:

library(rstan)

N <- 500

l1 <- 1.4; y1 <- rpois(N, l1)
l2 <- 0.8; y2 <- rpois(N, l2)
l3 <- 3; y3 <- rpois(N, l3)

l12 <- 1.2; y12 <- rpois(N, l12)
l13 <- 0.1; y13 <- rpois(N, l13)
l23 <- 0.5; y23 <- rpois(N, l23)

X <- matrix(0, N, 3)
X[,1] <- y1+y12+y13
X[,2] <- y2+y12+y23
X[,3] <- y3+y13+y23

model <- stan_model("MV_poisson.stan")
samples <- sampling(model, list(X=X,N=N), chains=1, iter=100)

Stan model code:

data {
  int N;
  int X[N, 3];
}

parameters {
  real<lower=0> lambda[6];
}

model {
  for(n in 1:N) {
      real lp[(1+min(X[n,1], X[n,2])) * (1+min(X[n,1], X[n,3])) * (1+min(X[n,2], X[n,3]))];
      int n_used;
      n_used <- 0;
      for(y12 in 0:min(X[n,1], X[n,2]))
       for(y13 in 0:min(X[n,1], X[n,3]))
          for(y23 in 0:min(X[n,2], X[n,3]))   
            if((y12+y13) <= X[n,1] && (y12+y23) <= X[n,2] && (y13+y23) <= X[n,3]) {
              n_used <- n_used + 1;
              lp[n_used] <- poisson_log(X[n,1]-y12-y13, lambda[1]) +
                          poisson_log(X[n,2]-y12-y23, lambda[2]) +
                          poisson_log(X[n,3]-y13-y23, lambda[3]) +
                          poisson_log(y12, lambda[4]) +
                          poisson_log(y13, lambda[5]) +
                          poisson_log(y23, lambda[6]);
            }
 
      increment_log_prob(log_sum_exp(lp[1:n_used]));
  }
 
  lambda ~ gamma(0.1, 0.1);
}

[andre.p...@googlemail.com]

Herra Huu,

thanks! Your solution works good. For a general solution, one could use

a constraint satisfaction problem (CSP) algorithm.

Andre

[Bob Carpenter]

Cool. Exactly what I meant by the marginalizing
being complicated --- I didn't have the heart to try
to code it!

- Bob

P.S. The X have to be integers so that the sums are integers.

> On Apr 21, 2016, at 4:41 AM, Herra Huu <heikki.uu...@gmail.com> wrote:
>
> As Aleksandar suspected, the current model doesn't correspond to MV Poisson distribution. In the current form it's just the product of marginal distributions, ie. defining joint density of form P(X1,X2,X3) = P(X1)P(X2)P(X3) and therefore missing the dependency structure.
>
> The correct joint density (for M=3) is given in page 36. Basically what you have to do is:
>

> take a sum over product of p(y1)p(y2)p(y3)p(y4)p(y5)p(y6), where each term is Poisson density, for all possible combinations of (y1,y2,y3,y12,y13,y23) such that:
> x1 = y1+y12+y13
> x2 = y2+y12+y23
> x3 = y3+y13+y23
>
> From there, it's easy to generalize for higher dimensions/different dependency structures, but it's also awfully computationally expensive if M/counts are high. Here is a simple example code for the M=3 case,