It’s not possible, but that’s by design rather than any mathematical restriction. The vectorised functions are designed to more efficiently compute the gradients and total likelihood for their inputs, so they’ll sample faster than a loop. The returned value is identical (within floating-point trickery) to taking the sum of a loop over each element.

So just to understand, do a vector of log-probabilities exist within the vectorized likelihood function at any point (but it’s returned as sum to the user), or does not exist?

If does exist, I would be interested in createing a C++ alternative where that vector of log-probabilities is returned.

do a vector of log-probabilities exist within the vectorized likelihood function at any point (but it’s returned as sum to the user), or does not exist

Generally no, but behaviour isn’t identical across all likelihoods. I believe the main issue is that the gradients are constructed in reference to returning the sum, and so both the return type and handling of gradients would need to change for all distributions which would be a pretty non-trivial change.

@bbbales2 Have I got that right? That for the vectorised distributions to return a vector of log-likelihoods rather than the single sum the current operands_and_partials approach would need a decent refactor?

The issue here is that returning each element instead of the sum creates a huge overhead internally. You would shoot yourself in the foot big time if you return things by element instead of the sum. It‘s probably even not such a performance difference to just do the loop in Stan as compared to working this out in C++.