There very well may be a clever way to do this faster, but my intuition is that this operation is going to be inherently slow because of the need to calculate the logp of a single data set associated with a single set of model parameters. But finding an efficient solution would be useful for the reasons you mention. Maybe @ricardoV94 or @junpenglao have some suggestions for speeding things up?