Comment 1. For instance, assume miny = 0. Now look at the term exp(-lambda1 - lambda2 - lambda3) in the equation, before the summation. In the code it says
real bipois_lpmf(int[] r ,real mu1, real mu2, real correlation_coeff) {
// r = argument to evaluate on
// https://en.wikipedia.org/wiki/Poisson_distribution#Bivariate_Poisson_distribution
// http://www2.stat-athens.aueb.gr/~karlis/multivariate%20Poisson%20models.pdf
real log_base_factor;
real log_sum_factor;
real log_factor;
int miny;
// Logarithm of the base factor in the multivariate Poisson distribution
// ie the term before the summation in the equation here
// https://en.wikipedia.org/wiki/Poisson_distribution#Bivariate_Poisson_distribution
log_base_factor = poisson_lpmf(r[1] | mu1) + poisson_lpmf(r[2] | mu2) -
correlation_coeff;
// Number of terms in the summation
miny = min(r[1], r[2]);
// This initial conditions works because the first term in the sum
// when k = 0 evaluates to 1, so
// log_sum_exp(0, log(term)) = log_sum ( 1 + term)
log_sum_factor = 0;
if(miny > 0) {
for(k in 1:miny) {
// The term is: choose(r_1, k) * choose(r_2, k) * k! * (corr / (mu1 * mu2))^k
// Here we compute the log of that term
log_factor = lchoose(r[1], k) + lchoose(r[2], k) +
lgamma(k + 1) +
k * (log(correlation_coeff) - log(mu1) - log(mu2));
// The log-sum-exp function is associative, since
// log_sum_exp(log_sum_exp(a, b), c) = log_sum_exp(a, b, c)
log_sum_factor = log_sum_exp(log_sum_factor, log_factor);
}
}
return log_base_factor + log_sum_factor;
}
Any comments on this? Am I onto something, or am I really dense or missing something obvious?
Hi!
Looking at bipois_lpmf, I cannot get it to match the Wikipedia equation for bivariate Poisson:
An alternative source is page 5 in this slide deck.
Comment 1. For instance, assume
miny = 0
. Now look at the termexp(-lambda1 - lambda2 - lambda3)
in the equation, before the summation. In the code it saysShould it not be
-mu3
instead of- exp(mu3)
since we're in log-space?Comment 2. Also, in the sum I don't understand where the terms in
come from. Seems wrong. When I take logs of the terms in the sum I get (see full code below):
Any comments on the above?
Here is my attempt at an implementation:
Any comments on this? Am I onto something, or am I really dense or missing something obvious?
Thanks in advance for your inputs.