# Loo and loglikelihood calculation for bivariate poisson

Given a standard model like y_i \sim \mathcal{Poisson}\left(\lambda\right), then the calculation of the log-likelihood function for using it with the loo package is straightforward:

vector[N] log_lik;
for (n in 1:N) {
log_lik[n] = poisson_lpmf(lambda);
}


My model is more complicated, since I have:

x_i \sim \mathcal{Poisson}\left( \lambda_1\right)\\ y_i \sim \mathcal{Poisson}\left( \lambda_2\right)

with x_i and y_i jointly distributed.
In this case, how can I write the log-likelihood calculation in the generated quantities of my model?
Thank you in advance!

This might be what you need

vector[N] log_lik;
for (n in 1:N) {
log_lik[n] = poisson_lpmf(x[n] | lambda1) + poisson_lpmf(y[n] | lambda2);
}


but I would need to see the model block to be certain

1 Like

Thanks @avehtari for the quick answer. I would like to compare the two models proposed by jake thompson to model soccer scores.
Here are the two models:

i) the bivariate Poisson model

model {
vector[num_games] lambda1;
vector[num_games] lambda2;
vector[num_games] lambda3;

// priors
raw_alpha ~ normal(0, 1);
raw_delta ~ normal(0, 1);
raw_rho ~ normal(0, 1);
mu ~ normal(0, 10);
eta ~ normal(0, 10);
sigma_a ~ normal(0, 10);
sigma_d ~ normal(0, 10);
sigma_r ~ normal(0, 10);

// likelihood
for (g in 1:num_games) {
lambda1[g] = exp(mu + (eta * homeg[g]) + alpha[home[g]] + delta[away[g]]);
lambda2[g] = exp(mu + alpha[away[g]] + delta[home[g]]);
lambda3[g] = exp(rho[home[g]] + rho[away[g]]);
}
h_goals ~ poisson(lambda1 + lambda3);
a_goals ~ poisson(lambda2 + lambda3);
}


ii) game random intercept model

model {
vector[num_games] lambda1;
vector[num_games] lambda2;

// priors
raw_alpha ~ normal(0, 1);                         // attacking random effects
raw_delta ~ normal(0, 1);                         // defending random effects
raw_gamma ~ normal(0, 1);                         // game random effects
mu ~ normal(0, 10);
eta ~ normal(0, 10);
sigma_a ~ normal(0, 10);
sigma_d ~ normal(0, 10);
sigma_g ~ normal(0, 10);

// likelihood
for (g in 1:num_games) {
lambda1[g] = exp(mu + (eta * homeg[g]) + alpha[home[g]] + delta[away[g]] + gamma[g]);
lambda2[g] = exp(mu + alpha[away[g]] + delta[home[g]] + gamma[g]);
}
h_goals ~ poisson(lambda1);
a_goals ~ poisson(lambda2);
}


As I expected, the Poisson terms are independent conditional on the parameters and thus the likelihood is product of those terms and thus log likelihood is sum of these terms. This

  h_goals ~ poisson(lambda1 + lambda3);
a_goals ~ poisson(lambda2 + lambda3);


can also be written as

  target += poisson_lpmf(h_goals | lambda1 + lambda3);
target += poisson_lpmf(a_goals | lambda2 + lambda3);


which could also be written as

  target += poisson_lpmf(h_goals | lambda1 + lambda3) + poisson_lpmf(a_goals | lambda2 + lambda3);


Given vector/array arguments, _lpmf and _lpdf functions sum together the log probabilities / densities, and to get the individual terms we need to use for loop in the generated quantities.

vector[N] log_lik;
for (n in 1:N) {
log_lik[n] = poisson_lpmf(h_goals[n] | lambda1 + lambda3) + poisson_lpmf(a_goals[n] | lambda2 + lambda3);
}


This will compute the joint likelihood given h_goals[n] and a_goals[n]. PSIS-LOO computes what would happen if nth log_lik term is removed form the target, and in this case it would correspond to leaving out both observations h_goals[n] and a_goals[n].

2 Likes

Thanks @avehtari for clarifying!