I do not think this will work. The thing is, I need to know which draw goes with which yrep and I guess you are only passing the yrep into the function, but not the iteration number as well which you would need.

The plots I am looking for are really simple, so maybe we can leave this out for bayesplot…unless you guys have clever ideas aside from pair-wise scatter plots.

What I have in mind to do for my given problem is to resample a subset of my fitted data (resampling data sets of the same size as my hold-out data set), then calculate T(y_subsampled, theta) and calculate the mean over the draws. I then apply this to a hod-out data set and compare the respective mean calculated on that hold out data set to my reference distribution. At least that is my plan… haven’t seen this somewhere else, but I think it makes sense for my application.