Bayesian approach to report simulation studies? I am running a simulation study where we want to estimate a proportion $p$. We are reporting the coverage of credible intervals with a uniform prior, and we are doing $500$ Monte Carlo simulations. We report the proportion of intervals that contain the true value of the parameter, and report a confidence interval for this proportion. I was wondering if I should also do a Bayesian analysis of the results from the simulation study or if the frequentist nature of the simulation study favours reporting confidence intervals.
Q. More in general, is it valid to use a Bayesian approach to report a simulation study? 
 A: The following is a simple proof of the (perhaps obvious) fact that a Bayesian credible interval (or set) always has average frequentist coverage equal to the nominal level (assuming a proper prior).  Some authors (e.g. Brown et al. 2001 and Fagerland, Lydersen & Laake 2017, section 1.4.2), in choosing between the large number of approximate frequentist confidence intervals, argue that such average coverage may be more relevant than ensuring that the coverage at least has the nominal level for all $\theta$.  If one accepts this, an average weighted by prior belief in different values of $\theta$ is arguably a sensible way to assess overall coverage (as done below).
By definition, a Bayesian $(1-\alpha)$-credible interval (or set) $C(x)$ (some function of the data $x$) includes $\theta$ with probability
$$
P(\theta \in C(x)|x)=1-\alpha. \tag{1}
$$
conditional on the data $x$.  Thinking of both $\theta$ and $x$ as random variables, it is also true, by the law of total probablity, that
$$
P(\theta \in C(x))=\int P(\theta \in C(x)|x)f(x)dx=1-\alpha \tag{2}
$$
where $f(x)$ is the marginal density of $x$.
For the frequentist coverage of the Bayesian credible interval, we instead condition on only the parameter $\theta$ such that $x$ and $C(x)$ are random.  This gives the frequentist coverage
$$
P(\theta \in C(x)|\theta) \tag{3}
$$
which in general will depend on $\theta$.  But the average frequentist coverage, weighted by the prior $\pi(\theta)$, is
$$
\int P(\theta \in C(x)|\theta) \pi(\theta) d\theta=P(\theta \in C(x)) =1-\alpha, \tag{4}
$$
by the law of total probability and (2), and so always perfectly matches the nominal level.  
The following plot and R-code illustrates (3) and (4) for a binomial model with a uniform prior on $p$ and $n=20$:

R code
m <- 1e+5
p <- seq(0,1,len=m)
coverage <- numeric(m)
alpha <- .05
a <- 1
b <- 1
n <- 20
x <- 0:n # possible outcomes
lower <- qbeta(alpha/2, a+x, b+n-x) # lower and upper credible limits
upper <- qbeta(1-alpha/2, a+x, b+n-x) # for each value of x
for (i in 1:length(p)) {
  px <- dbinom(x,size=n,prob=p[i]) # probability of possible outcomes given p
  coverage[i] <- sum(px*ifelse(lower<=p[i] & p[i]<=upper, 1,0)) # coverage
}
plot(p,coverage,type="l",ylim=c(.85,1))
abline(h=1-alpha, lty=3)

> mean(coverage)
[1] 0.9499928

A: This is probably a bit opinion based, but I would normally expect that many of the usual arguments for using Bayesian methods may not apply. E.g. bringing in extra information to compensate for low sample size (usually one can just simulate more) and wanting to use prior information (firstly, exact prior information may be hard to come but and it may be less controversial to just simulate more) would not seem like strong arguments. 
On the other hand, using a Bayesian method to get an easy to compute approximation to (exact) frequentist confidence intervals using percentiles of the beta distribution may be a reasonable idea - as long as you have enough simulations that your exact choice of prior will be uncontroversial. Putting too strong a prior on a desirable result (e.g. using a Beta(95,5) prior for the coverage probability of 95% credible intervals) might end up being hard to defend, while I would guess most people would not have much of a problem using vague priors such as e.g. Beta(1/3, 1/3), Beta(0.5, 0.5) or the uniform Beta(1, 1) prior, if you have a decent number of simulations.
The obvious frequentist alternative is something like an exact Clopper-Pearson confidence interval. I have never had anyone questioning (e.g. on the basis of methodological incoherence) me using one those for simulation studies, even when I was evaluating Bayesian methods. In a way, you are after all evaluating the frequentist properties of a method (repeat sampling performance). Around that I dimly recall some quote (perhaps by Rubin?) that you use the Bayesian approach to construct methods and evaluate these properties in a frequentist way (and that they usually turn out to have good frequentist properties).
On the specifics of your project, 500 simulations seems like a pretty low number for the coverage probability of credible intervals (unless they are e.g. 50% intervals, certainly quite low for 95% ones). Even if you observe exactly 95% of simulations with the CI covering your true value, your credible interval for the credible interval coverage under a Beta(0.5, 0.5) prior would be from 92.8% to 96.7% (2500 simulations would get you to 94.1% to 95.8% and 10,000 to 94.6 to 95.4%). I guess if that level of uncertainty is fine in your case that may not be an issue.
