Loo and loglikelihood calculation for bivariate poisson

Modeling
#1

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!

#2

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
#3

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);
}
#4

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
#5

Thanks @avehtari for clarifying!