laplace_marginal: inconsistency between doc and math for built-in bernoulli_logit likelihood

[charlesm93]

This concerns laplace_marginal_bernoulli_logit and related functions.

Our current doc (https://mc-stan.org/docs/functions-reference/embedded_laplace.html) indicates that users should pass y, a vector of binary observations and y_index, an array which indicates to which group each y belongs to (i.e., which element of theta y_i corresponds to).

However the C++ code in stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp uses a different implementation which requires users to pass the sum of all observations, $\sum_{i \in g(i)} y_i$ for each group and the number of elements in each group, $\sum_{i \in g(i)} 1$, which are sufficient statistics for the likelihood.

The implementation in C++ corresponds to what I had originally built for the prototype embedded Laplace approximation. However, we decided to move to the signature in the doc and we've done this correctly for the Poisson and NegativeBinomial likelhoods, but it looks like Bernoulli_logit slipped through the cracks....

This should be a straightforward fix.

[SteveBronder]

I think this is fine and reasonable. Though for the signature where the user has to pass in just the vector of binary observations y and the array with indicates which group each y belongs to there will be extra work to calculate the counts per group each time. For example, this is our poisson likelihood where you can see we do the sum internally.

struct poisson_log_likelihood {
  template <typename Theta, typename YVec, typename Mean,
            require_eigen_vector_t<Theta>* = nullptr>
  inline auto operator()(const Theta& theta, const YVec& y,
                         const std::vector<int>& y_index, Mean&& mean,
                         std::ostream* /*pstream*/) const {
    Eigen::VectorXd counts_per_group = Eigen::VectorXd::Zero(theta.size());
    Eigen::VectorXd n_per_group = Eigen::VectorXd::Zero(theta.size());

    for (int i = 0; i < theta.size(); i++) {
      counts_per_group(y_index[i] - 1) += y[i];
      n_per_group(y_index[i] - 1) += 1;
    }

    auto theta_offset = to_ref(add(theta, mean));

    return -sum(lgamma(add(counts_per_group, 1)))
           + dot_product(theta_offset, counts_per_group)
           - dot_product(n_per_group, exp(theta_offset));
  }
};

Since y and y_index have to be std:vector<int> types here counts_per_group and n_per_group could have been calculated by the user once in the transformed data block which would save us this computation for each likelihood evaluation.

What is the reasoning behind preferring the y and y_index version over the pre summarized version?

We could have both signatures

laplace_marginal_poisson_log_lpmf(y, y_index,...)` 

and

laplace_marginal_poisson_log_summary_lpmf(counts_per_group, n_per_group, ...)

[charlesm93]

Yes, you're correct that this is more calculation and that's why my initial implementation of the Poisson and Bernoulli likelihoods used the sum per group and size of group as inputs.

While less efficient, taking y and y_index is more user-friendly and I think that's why we made the switch for the Poisson likelihood. I'm not too keen on adding another wrapper and I haven't found the efficiency gain to be important. I think if this really matters to the user, they can specify their own likelihood with laplace_marginal.

Maybe @avehtari can weight in.

[charlesm93]

p.s. I saw you self-assigned the issue. I'm also happy to handle this one, but up to you.