Hello Stan users. I am new to Stan (I started using Stan today) and I was wondering what is the correct way of setting up a compound truncated distribution. I am working with count data (the list consists of non-zero positive integers) and I want to fit a zero-truncated Poisson Lognormal (PLN) distribution to it. I know that the Poisson Lognormal distribution is just a Poisson distribution with a rate parameter that is log-normally distributed.
data{
int<lower=0> N;
int y[N];
}
parameters {
real mu;
real<lower=0> sigma;
vector[N] log_lambda;
}
model {
sigma ~ cauchy(0, 5);
mu ~ normal(0, 10);
log_lambda ~ normal(mu, sigma);
for (n in 1:N)
target += poisson_log_lpmf(y[n] | log_lambda[n]);
}
How do I properly modify it to account for the truncation? I am still having a hard time understanding what the target is and why does this code work. Thank you very much!
First of all, I assume that your data would not include zeroes, correct? This is not reflected in the declaration of data y.
To model the data as coming from a zero-truncated Poisson, you essentially just have to renormalise the posterior to have area under the curve of 1 after truncation of zero. On the probability scale this is achieved by dividing by one minus the probability of observing a zero count under the non-truncated distribution. Stan works on the log probability scale, so that means you need to subtract the logarithm of 1 minus the probability of observing a zero, given the Poisson rate.
Now, as you are modeling a separate (log) rate for each count observation, you will have to subtract log(1 - poisson_log_lpmf(0 | log_lambda[n])); from the log posterior (i.e. “target”) for each of the N observations. Note that there is a numerically stable implementation of log(1-x) in Stan, so look that up and use that one rather than my pseudo code above.
Hello. Thank you for your quick reply. Yes I am working with a vector that has counts greater than zero. I have two questions about this. For the truncation term, wouldn’t I have to exponentiate the distribution evaluated at zero? That is log(1 - exp(poisson_log_lpmf(0|log_lambda[n])) or simply log1m_exp(poisson_log_lpmf(0|log_lambda[n])) for the numerically stable version.
The other question was… for any compound model (like the PLN that must be integrated over all possible values of the rate parameter lambda), should I always account for this integration by setting a different lambda for each data point like it is done in the sample code? I don’t exactly understand why this must be true.
You’re absolutely right. Typo on my part! Note that log1m(poisson_log_pmf(0|log_lambda[n])also suffices. (so just pmf and not lpmf at the end of the poisson statement).
Yes you need to always integrate over the mixture component in compound distributions. In this particular case there is also only one way to do it as there is no analytical form for the integration over lambda. In your Stan model, you therefore have to specify a separate lambda for each count and integrate over all potential values via the Monte Carlo sampling, which we for Stan is fortunately straightforward as lambda is a continuous variable.
For some other compound distributions, there is an analytical solutions to integrating over the mixture component, like for the gamma-Poisson compound distribution, where lambda is assumed to be gamma-distributed, and the integral over lambda yields a negative binomial distribution with mean and shape equal to the mean and shape of the gamma mixture component - see eg wikipedia for some more explanation).